From 392d1d93d87eb18d9fd721ccfc848db6d04c5e3d Mon Sep 17 00:00:00 2001 From: xie7 Date: Wed, 2 Oct 2024 14:01:15 -0700 Subject: [PATCH] Subject (<= 72 Characters) Added orographic drag toolkit for topo generation. *********1*********2*********3*********4*********5*********6*********7** Longer commit message body describing the commit. Can contain lists as follows: * Item 1 * Item 2 * Item 3 A good commit message should be written like an email, a subject followed by a blank line, followed by a more descriptive body. Can also contain a tag at the bottom describing what type of commit this is. [BFB] - Bit-For-Bit [FCC] - Flag Climate Changing [Non-BFB] - Non Bit-For-Bit [CC] - Climate Changing [NML] - Namelist Changing See confluence for a more detailed description about these tags. --- .../tools/orographic_drag_toolkit/Makefile | 106 + .../eam/tools/orographic_drag_toolkit/README | 18 + .../Tempest-remap_generation.sh | 13 + .../cube_to_target.F90 | 2680 +++++++++++++++++ .../tools/orographic_drag_toolkit/make.ncl | 21 + .../orographic_drag_toolkit/ogwd_sub.F90 | 908 ++++++ .../orographic_drag_toolkit/reconstruct.F90 | 2675 ++++++++++++++++ .../tools/orographic_drag_toolkit/remap.F90 | 1564 ++++++++++ .../eam/tools/orographic_drag_toolkit/run.sh | 6 + .../orographic_drag_toolkit/shr_kind_mod.F90 | 20 + .../orographic_drag_toolkit/transform.F90 | 351 +++ 11 files changed, 8362 insertions(+) create mode 100755 components/eam/tools/orographic_drag_toolkit/Makefile create mode 100755 components/eam/tools/orographic_drag_toolkit/README create mode 100755 components/eam/tools/orographic_drag_toolkit/Tempest-remap_generation.sh create mode 100755 components/eam/tools/orographic_drag_toolkit/cube_to_target.F90 create mode 100755 components/eam/tools/orographic_drag_toolkit/make.ncl create mode 100755 components/eam/tools/orographic_drag_toolkit/ogwd_sub.F90 create mode 100755 components/eam/tools/orographic_drag_toolkit/reconstruct.F90 create mode 100755 components/eam/tools/orographic_drag_toolkit/remap.F90 create mode 100755 components/eam/tools/orographic_drag_toolkit/run.sh create mode 100755 components/eam/tools/orographic_drag_toolkit/shr_kind_mod.F90 create mode 100755 components/eam/tools/orographic_drag_toolkit/transform.F90 diff --git a/components/eam/tools/orographic_drag_toolkit/Makefile b/components/eam/tools/orographic_drag_toolkit/Makefile new file mode 100755 index 000000000000..ec236185cf67 --- /dev/null +++ b/components/eam/tools/orographic_drag_toolkit/Makefile @@ -0,0 +1,106 @@ +EXEDIR = . +EXENAME = cube_to_target +RM = rm + +.SUFFIXES: +.SUFFIXES: .F90 .o + +FC = ifort +DEBUG = FALSE + + +# Check for the NetCDF library and include directories +ifeq ($(LIB_NETCDF),$(null)) +LIB_NETCDF := /gpfs/fs1/soft/chrysalis/spack/opt/spack/linux-centos8-x86_64/intel-20.0.4/netcdf-fortran-4.4.4-rdxohvp/lib +#/global/common/software/nersc/pm-2023q1/spack-stacks-1/views/climate-utils/lib +#/public/software/mathlib/netcdf/4.3.2/intel/lib +endif + +ifeq ($(INC_NETCDF),$(null)) +INC_NETCDF := /gpfs/fs1/soft/chrysalis/spack/opt/spack/linux-centos8-x86_64/intel-20.0.4/netcdf-fortran-4.4.4-rdxohvp/include +#/global/common/software/nersc/pm-2023q1/spack-stacks-1/views/climate-utils/include +#/public/software/mathlib/netcdf/4.3.2/intel/include +endif + +# Determine platform +UNAMES := $(shell uname -s) +UNAMEM := $(findstring CRAY,$(shell uname -m)) + +#------------------------------------------------------------------------ +# LF95 +#------------------------------------------------------------------------ +# +# setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:/usr/local/netcdf-4.1.3-gcc-4.4.4-13-lf9581/lib +# +ifeq ($(FC),lf95) +# +# Tramhill +# + INC_NETCDF :=/usr/local/netcdf-4.1.3-gcc-4.4.4-13-lf9581/include + LIB_NETCDF :=/usr/local/netcdf-4.1.3-gcc-4.4.4-13-lf9581/lib + + LDFLAGS = -L$(LIB_NETCDF) -lnetcdf -lnetcdff -lcurl -lhdf5 -lhdf5_hl -mcmodel=medium + FFLAGS := -c --trace --trap --wide -CcdRR8 -I$(INC_NETCDF) + ifeq ($(DEBUG),TRUE) +# FFLAGS += --chk aesu -Cpp --trace + FFLAGS += -g --chk a,e,s,u --pca + else + FFLAGS += -O + endif + +endif + + + +.F90.o: + $(FC) $(FFLAGS) $< + + +#------------------------------------------------------------------------ +# AIX +# #------------------------------------------------------------------------ +# + #ifeq ($(UNAMES),AIX) + FC = ifort #xlf90 + #FFLAGS = -c -I$(INC_NETCDF) -I/BIGDATA1/iapcas_mhzhang_xiejinbo/topo_tool/cube_to_target/functional/ -convert big_endian + + FFLAGS = -c -I$(INC_NETCDF) -convert big_endian -traceback + #FFLAGS := -c -I$(INC_NETCDF) -no-prec-div -traceback -convert big_endian -fp-model source -assume byterecl -ftz -m64 -mcmodel=large -safe-cray-ptr + LDFLAGS = -L$(LIB_NETCDF) -lnetcdff + #LDFLAGS = -L$(LIB_NETCDF) -lnetcdf -lnetcdff -m64 -static-intel + .F90.o: + $(FC) $(FFLAGS) -qsuffix=f=F90 $< +# #endif + + +.F90.o: + $(FC) $(FFLAGS) $< + + + + + + + + + + +#------------------------------------------------------------------------ +# Default rules and macros +#------------------------------------------------------------------------ + +#OBJS := reconstruct.o remap.o cube_to_target.o shr_kind_mod.o +OBJS := reconstruct.o remap.o shr_kind_mod.o transform.o sub_xjb.o cube_to_target.o +#OBJS := reconstruct.o remap.o cube_to_target.o sub.o shr_kind_mod.o +#sub.o shr_kind_mod.o + +$(EXEDIR)/$(EXENAME): $(OBJS) + $(FC) -o $@ $(OBJS) -I$(INC_NETCDF) $(LDFLAGS) + +clean: + $(RM) -f $(OBJS) *.mod $(EXEDIR)/$(EXENAME) + +cube_to_target.o: shr_kind_mod.o remap.o reconstruct.o sub_xjb.o transform.o +remap.o: +reconstruct.o: remap.o +#reconstruct.o : shr_kind_mod.o diff --git a/components/eam/tools/orographic_drag_toolkit/README b/components/eam/tools/orographic_drag_toolkit/README new file mode 100755 index 000000000000..1675a91d5e76 --- /dev/null +++ b/components/eam/tools/orographic_drag_toolkit/README @@ -0,0 +1,18 @@ +cube_to_target performs rigourous remapping of topo variables from cubed-sphere grid to +any target grid. In the process SGH is computed. + +Input files: + +1. USGS-topo-cube.nc (may be found here $CESMDATA/inputdata/atm/cam/hrtopo/USGS-topo-cube3000.nc) + + This is the topo data on a cubed-sphere (default is 3km cubed-sphere grid) + +2. target.nc (e.g., $CESMDATA/inputdata/atm/cam/grid-description/se/ne30np4_091226_pentagons.nc) + + This is a SCRIP/ESMF grid descriptor file for the target grid + +3. phis-smooth.nc + + (optional) The user may provide a smoothed PHIS field. The software then recomputes SGH to + account for the smoothing in the sub-grid-scale. + diff --git a/components/eam/tools/orographic_drag_toolkit/Tempest-remap_generation.sh b/components/eam/tools/orographic_drag_toolkit/Tempest-remap_generation.sh new file mode 100755 index 000000000000..e9bb8470393d --- /dev/null +++ b/components/eam/tools/orographic_drag_toolkit/Tempest-remap_generation.sh @@ -0,0 +1,13 @@ + + +source /lcrc/soft/climate/e3sm-unified/load_latest_e3sm_unified_chrysalis.sh +tempest_root=~/.conda/envs/jinbo +# Generate the element mesh. +${tempest_root}/bin/GenerateCSMesh --alt --res 30 --file topo2/ne30.g +# Generate the target physgrid mesh. +${tempest_root}/bin/GenerateVolumetricMesh --in topo2/ne30.g --out topo2/ne30pg2.g --np 2 --uniform +# Generate a high-res target physgrid mesh for cube_to_target. +${tempest_root}/bin/GenerateVolumetricMesh --in topo2/ne30.g --out topo2/ne30pg4.g --np 4 --uniform +# Generate SCRIP files for cube_to_target. +${tempest_root}/bin/ConvertMeshToSCRIP --in topo2/ne30pg4.g --out topo2/ne30pg4_scrip.nc +${tempest_root}/bin/ConvertMeshToSCRIP --in topo2/ne30pg2.g --out topo2/ne30pg2_scrip.nc diff --git a/components/eam/tools/orographic_drag_toolkit/cube_to_target.F90 b/components/eam/tools/orographic_drag_toolkit/cube_to_target.F90 new file mode 100755 index 000000000000..0bb64273b088 --- /dev/null +++ b/components/eam/tools/orographic_drag_toolkit/cube_to_target.F90 @@ -0,0 +1,2680 @@ +! +! DATE CODED: Nov 7, 2011 to Oct 15, 2012 +! DESCRIPTION: Remap topo data from cubed-sphere grid to target grid using rigorous remapping +! (Lauritzen, Nair and Ullrich, 2010, J. Comput. Phys.) +! +! Author: Peter Hjort Lauritzen (pel@ucar.edu), AMP/CGD/NESL/NCAR +! +program convterr + use shr_kind_mod, only: r8 => shr_kind_r8 + use reconstruct + use ogwd_sub + implicit none +# include + + !************************************** + ! + ! USER SETTINGS BELOW + ! + !************************************** + ! + ! + ! if smoothed PHIS is available SGH needs to be recomputed to account for the sub-grid-scale + ! variability introduced by the smoothing + ! +logical :: lsmooth_terr = .FALSE. +!logical :: lsmooth_terr = .TRUE. + ! + ! PHIS is smoothed by other software/dynamical core + ! + logical :: lexternal_smooth_terr = .FALSE. ! lexternal_smooth_terr = .FALSE. is NOT supported currently +!logical :: lexternal_smooth_terr = .TRUE. + ! + ! set PHIS=0.0 if LANDFRAC<0.01 + ! + logical :: lzero_out_ocean_point_phis = .TRUE.!.FALSE. +!logical :: lzero_out_ocean_point_phis = .FALSE. + ! + ! For internal smoothing (experimental at this point) + ! =================================================== + ! + ! if smoothing is internal (lexternal_smooth_terr=.FALSE.) choose coarsening factor + ! + ! recommendation: 2*(target resolution)/(0.03 degree) + ! + ! factor must be an even integer + ! + integer, parameter :: factor = 60 !coarse grid = 2.25 degrees + integer, parameter :: norder = 2 + integer, parameter :: nmono = 0 + integer, parameter :: npd = 1 + ! + !********************************************************************** + ! + ! END OF USER SETTINS BELOW + ! (do not edit beyond this point unless you know what you are doing!) + ! + !********************************************************************** + ! + integer :: im, jm, ncoarse + integer :: ncube !dimension of cubed-sphere grid + + real(r8), allocatable, dimension(:) :: landm_coslat, landfrac, terr, sgh30 + real(r8), allocatable, dimension(:) :: terr_coarse !for internal smoothing + + integer :: alloc_error,dealloc_error + integer :: i,j,n,k,index + integer*2, allocatable, dimension(:,:) :: iterr ! terrain data for 30-sec tile + integer ncid,status, dimlatid,dimlonid, landid, topoid ! for netCDF USGS data file + integer :: srcid,dstid, jm_dbg ! for netCDF weight file + integer, dimension(2) :: src_grid_dims ! for netCDF weight file + + integer :: dimid + + logical :: ldbg + real(r8), allocatable, dimension(:) :: lon , lat + real(r8), allocatable, dimension(:) :: lon_landm , lat_landm + real(r8), allocatable, dimension(:) :: area + integer :: im_landm, jm_landm + integer :: lonid, latid, phisid + ! + ! constants + ! + REAL (r8), PARAMETER :: pi = 3.14159265358979323846264338327 + REAL (r8), PARAMETER :: piq = 0.25*pi + REAL (r8), PARAMETER :: pih = 0.50*pi + REAL (r8), PARAMETER :: deg2rad = pi/180.0 + + real(r8) :: wt,dlat + integer :: ipanel,icube,jcube + real(r8), allocatable, dimension(:,:,:) :: weight,terr_cube,landfrac_cube,sgh30_cube + real(r8), allocatable, dimension(:,:,:) :: landm_coslat_cube + integer, allocatable, dimension(:,:) :: idx,idy,idp + integer :: npatch, isub,jsub, itmp, iplm1,jmin,jmax + real(r8) :: sum,dx,scale,dmax,arad,jof,term,s1,c1,clon,iof,dy,s2,c2,dist + ! + ! for linear interpolation + ! + real(r8) :: lambda,theta,wx,wy,offset + integer :: ilon,ilat,ip1,jp1 + ! + ! variable for regridding + ! + integer :: src_grid_dim ! for netCDF weight file + integer :: n_a,n_b,n_s,n_aid,n_bid,n_sid + integer :: count + real(r8), allocatable, dimension(:) :: landfrac_target, terr_target, sgh30_target, sgh_target +!====Jinbo Xie===== + real(r8), allocatable, dimension(:) :: oc_target + real(r8), allocatable, dimension(:,:) :: oa_target,ol_target + real(r8) :: terr_if + real(r8), allocatable, dimension(:) :: lat_terr,lon_terr + integer :: nvar_dirOA,nvar_dirOL + integer,allocatable,dimension(:) :: indexb !max indice dimension + real(r8),allocatable,dimension(:,:,:) :: terrout + real(r8),allocatable,dimension(:,:) :: dxy +!====Jinbo Xie===== + + real(r8), allocatable, dimension(:) :: landm_coslat_target, area_target + ! + ! this is only used if target grid is a lat-lon grid + ! + integer , parameter :: im_target = 360 , jm_target = 180 + ! + ! this is only used if target grid is not a lat-lon grid + ! + real(r8), allocatable, dimension(:) :: lon_target, lat_target + ! + ! new + ! + integer :: ntarget, ntarget_id, ncorner, ncorner_id, nrank, nrank_id + integer :: ntarget_smooth + real(r8), allocatable, dimension(:,:):: target_corner_lon, target_corner_lat + real(r8), allocatable, dimension(:) :: target_center_lon, target_center_lat, target_area +!==========Jinbo Xie============ +real(r8), allocatable, dimension(:,:):: target_corner_lon_deg,target_corner_lat_deg +!==========Jinbo Xie============ + integer :: ii,ip,jx,jy,jp + real(r8), dimension(:), allocatable :: xcell, ycell, xgno, ygno + real(r8), dimension(:), allocatable :: gauss_weights,abscissae + integer, parameter :: ngauss = 3 + integer :: jmax_segments,jall + real(r8) :: tmp + + real(r8), allocatable, dimension(:,:) :: weights_all + integer , allocatable, dimension(:,:) :: weights_eul_index_all + integer , allocatable, dimension(:) :: weights_lgr_index_all + integer :: ix,iy + ! + ! volume of topography + ! + real(r8) :: vol_target, vol_target_un, area_target_total,vol_source,vol_tmp + integer :: nlon,nlon_smooth,nlat,nlat_smooth + logical :: ltarget_latlon,lpole + real(r8), allocatable, dimension(:,:) :: terr_smooth + ! + ! for internal filtering + ! + real(r8), allocatable, dimension(:,:) :: weights_all_coarse + integer , allocatable, dimension(:,:) :: weights_eul_index_all_coarse + integer , allocatable, dimension(:) :: weights_lgr_index_all_coarse + real(r8), allocatable, dimension(:) :: area_target_coarse + real(r8), allocatable, dimension(:,:) :: da_coarse,da + real(r8), allocatable, dimension(:,:) :: recons,centroids + integer :: nreconstruction + + integer :: jmax_segments_coarse,jall_coarse,ncube_coarse + real(r8) :: all_weights + character(len=512) :: target_grid_file + character(len=512) :: input_topography_file + character(len=512) :: output_topography_file + character(len=512) :: smoothed_topography_file +!=======Jinbo Xie array======= +real(r8) :: xxt,yyt,zzt +!real(r8),allocatable,dimension(:) :: xbar,ybar,zbar + +real(r8),dimension(32768) :: xhds,yhds,zhds,hds,xbar,ybar,zbar,lon_bar,lat_bar +real(r8) :: rad,xx2,yy2,zz2,ix2,iy2,ip2 +real(r8) :: lonii,latii +character*20 :: indice + +!=======Jinbo Xie array======= +nvar_dirOA=2+1!4 !2+1!4!36 +nvar_dirOL=180 +!=======Jinbo Xie array======= + + ! + ! turn extra debugging on/off + ! + ldbg = .FALSE. + + nreconstruction = 1 + ! + call parse_arguments(target_grid_file , input_topography_file , & + output_topography_file, smoothed_topography_file, & + lsmooth_terr ) + ! + !********************************************************* + ! + ! read in target grid + ! + !********************************************************* + ! + status = nf_open(trim(target_grid_file), 0, ncid) + IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR(STATUS) + + status = NF_INQ_DIMID(ncid, 'grid_size', ntarget_id) + status = NF_INQ_DIMLEN(ncid, ntarget_id, ntarget) + WRITE(*,*) "dimension of target grid: ntarget=",ntarget + + status = NF_INQ_DIMID(ncid, 'grid_corners', ncorner_id) + status = NF_INQ_DIMLEN(ncid, ncorner_id, ncorner) + WRITE(*,*) "maximum number of corners: ncorner=",ncorner + + status = NF_INQ_DIMID(ncid, 'grid_rank', nrank_id);status = NF_INQ_DIMLEN(ncid, nrank_id, nrank) + WRITE(*,*) "grid rank: nrank=",nrank + IF (nrank==2) THEN + WRITE(*,*) "target grid is a lat-lon grid" + ltarget_latlon = .TRUE. + status = NF_INQ_DIMID(ncid, 'nlon', ntarget_id) + status = NF_INQ_DIMLEN(ncid, ntarget_id, nlon) + status = NF_INQ_DIMID(ncid, 'nlat', ntarget_id) + status = NF_INQ_DIMLEN(ncid, ntarget_id, nlat) + status = NF_INQ_DIMID(ncid, 'lpole', ntarget_id) + status = NF_INQ_DIMLEN(ncid, ntarget_id, lpole) +!=====Jinbo Xie====== +!nlon=256!720!256 +!nlat=128!361!128 +!nvar_dirOA=2+1!4 !2+1!4!36 +!nvar_dirOL=180 !4!8 !4 !720 !180 !4 !180 !4!180 !720 !180 !4 !180 !720!180 !4 !180 !4 !2*11520!4!720!5760!2880!1440!720 !360!180 !4 !180!360!4 !18!36!4 + +!=====Jinbo Xie====== + WRITE(*,*) "nlon=",nlon,"nlat=",nlat + IF (lpole) THEN + WRITE(*,*) "center of most Northern grid cell is lat=90; similarly for South pole" + ELSE + WRITE(*,*) "center of most Northern grid cell is NOT lat=90; similarly for South pole" + END IF + ELSE IF (nrank==1) THEN + ltarget_latlon = .FALSE. + ELSE + WRITE(*,*) "nrank out of range",nrank + STOP + ENDIF + + allocate ( target_corner_lon(ncorner,ntarget),stat=alloc_error) + allocate ( target_corner_lat(ncorner,ntarget),stat=alloc_error) + !====Jinbo Xie====== + allocate ( target_corner_lon_deg(ncorner,ntarget),stat=alloc_error) + allocate ( target_corner_lat_deg(ncorner,ntarget),stat=alloc_error) + !====Jinbo Xie====== + status = NF_INQ_VARID(ncid, 'grid_corner_lon', lonid) + status = NF_GET_VAR_DOUBLE(ncid, lonid,target_corner_lon) +!Jinbo Xie + target_corner_lon_deg=target_corner_lon +!Jinbo Xie + IF (maxval(target_corner_lon)>10.0) target_corner_lon = deg2rad*target_corner_lon + + status = NF_INQ_VARID(ncid, 'grid_corner_lat', latid) + status = NF_GET_VAR_DOUBLE(ncid, latid,target_corner_lat) +!Jinbo Xie + target_corner_lat_deg=target_corner_lat +!Jinbo Xie + IF (maxval(target_corner_lat)>10.0) target_corner_lat = deg2rad*target_corner_lat + ! + ! for writing remapped data on file at the end of the program + ! + allocate ( target_center_lon(ntarget),stat=alloc_error) + allocate ( target_center_lat(ntarget),stat=alloc_error) + allocate ( target_area (ntarget),stat=alloc_error)!dbg + + status = NF_INQ_VARID(ncid, 'grid_center_lon', lonid) + status = NF_GET_VAR_DOUBLE(ncid, lonid,target_center_lon) + + status = NF_INQ_VARID(ncid, 'grid_center_lat', latid) + status = NF_GET_VAR_DOUBLE(ncid, latid,target_center_lat) + + status = NF_INQ_VARID(ncid, 'grid_area', latid) + status = NF_GET_VAR_DOUBLE(ncid, latid,target_area) + + status = nf_close (ncid) + if (status .ne. NF_NOERR) call handle_err(status) + ! + !**************************************************** + ! + ! get dimension of cubed-sphere grid + ! + !**************************************************** + ! + WRITE(*,*) "get dimension of cubed-sphere data from file" + !status = nf_open('USGS-topo-cube3000.nc', 0, ncid) + status = nf_open(trim(input_topography_file), 0, ncid) + IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR(STATUS) + status = NF_INQ_DIMID(ncid, 'grid_size', dimid) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + status = NF_INQ_DIMLEN(ncid, dimid, n) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + + ncube = INT(SQRT(DBLE(n/6))) + WRITE(*,*) "cubed-sphere dimension: ncube = ",ncube + WRITE(*,*) "average grid-spacing at the Equator (degrees):" ,90.0/ncube + + status = nf_close (ncid) + if (status .ne. NF_NOERR) call handle_err(status) + ! + !**************************************************** + ! + ! compute weights for remapping + ! + !**************************************************** + ! + jall = ncube*ncube*12*10 !anticipated number of weights (cab be tweaked) + jmax_segments = 100000 !can be tweaked + + allocate (weights_all(jall,nreconstruction),stat=alloc_error ) + allocate (weights_eul_index_all(jall,3),stat=alloc_error ) + allocate (weights_lgr_index_all(jall),stat=alloc_error ) + +!!======Jinbo Xie==== +!!Jinbo Xie debug +!#if 0 + CALL overlap_weights(weights_lgr_index_all,weights_eul_index_all,weights_all,& + jall,ncube,ngauss,ntarget,ncorner,jmax_segments,target_corner_lon,target_corner_lat,nreconstruction) +!#endif +!weights_all=0.01 +!!Jinbo Xie debug +!!======Jinbo Xie====== + ! + !**************************************************** + ! + ! read cubed-sphere 3km data + ! + !**************************************************** + ! + WRITE(*,*) "read cubed-sphere 3km data from file" + status = nf_open('USGS-topo-cube3000.nc', 0, ncid) + IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR(STATUS) + + status = NF_INQ_DIMID(ncid, 'grid_size', dimid) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + status = NF_INQ_DIMLEN(ncid, dimid, n) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + + ncube = INT(SQRT(DBLE(n/6))) + WRITE(*,*) "cubed-sphere dimension, ncube: ",ncube + + allocate ( landm_coslat(n),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for landfrac' + stop + end if + + status = NF_INQ_VARID(ncid, 'LANDM_COSLAT', landid) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + + status = NF_GET_VAR_DOUBLE(ncid, landid,landm_coslat) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + WRITE(*,*) "min/max of landm_coslat",MINVAL(landm_coslat),MAXVAL(landm_coslat) + ! + ! read LANDFRAC + ! + allocate ( landfrac(n),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for landfrac' + stop + end if + + status = NF_INQ_VARID(ncid, 'LANDFRAC', landid) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + + status = NF_GET_VAR_DOUBLE(ncid, landid,landfrac) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + WRITE(*,*) "min/max of landfrac",MINVAL(landfrac),MAXVAL(landfrac) + ! + ! read terr + ! + allocate ( terr(n),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for landfrac' + stop + end if + + status = NF_INQ_VARID(ncid, 'terr', landid) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + + status = NF_GET_VAR_DOUBLE(ncid, landid,terr) + + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + WRITE(*,*) "min/max of terr",MINVAL(terr),MAXVAL(terr) +!===Jinbo Xie read in lat lon==== + allocate ( lat_terr(n),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for lat_terr' + stop + end if + status = NF_INQ_VARID(ncid, 'lat', landid) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + status = NF_GET_VAR_DOUBLE(ncid, landid,lat_terr) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + WRITE(*,*) "min/max of lat",MINVAL(lat_terr),MAXVAL(lat_terr) + + allocate ( lon_terr(n),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for lon_terr' + stop + end if + status = NF_INQ_VARID(ncid, 'lon', landid) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + + status = NF_GET_VAR_DOUBLE(ncid, landid,lon_terr) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + WRITE(*,*) "min/max of lon",MINVAL(lon_terr),MAXVAL(lon_terr) +!===Jinbo Xie read in lat lon==== + ! + ! + ! + allocate ( sgh30(n),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for landfrac' + stop + end if + + status = NF_INQ_VARID(ncid, 'SGH30', landid) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + + status = NF_GET_VAR_DOUBLE(ncid, landid,sgh30) + IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) + WRITE(*,*) "min/max of sgh30",MINVAL(sgh30),MAXVAL(sgh30) + + print *,"close file" + status = nf_close (ncid) + if (status .ne. NF_NOERR) call handle_err(status) + + WRITE(*,*) 'done reading in LANDM_COSLAT data from netCDF file' + ! + !********************************************************* + ! + ! do actual remapping + ! + !********************************************************* + ! + allocate (terr_target(ntarget),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for terr_target' + stop + end if + allocate (landfrac_target(ntarget),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for landfrac_target' + stop + end if + allocate (landm_coslat_target(ntarget),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for landfrac_target' + stop + end if + allocate (sgh30_target(ntarget),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for sgh30_target' + stop + end if + allocate (area_target(ntarget),stat=alloc_error ) + terr_target = 0.0 + landfrac_target = 0.0 + sgh30_target = 0.0 + landm_coslat_target = 0.0 + area_target = 0.0 + + tmp = 0.0 + do count=1,jall + i = weights_lgr_index_all(count) + wt = weights_all(count,1) + area_target (i) = area_target(i) + wt + end do + + +!#if 0 +!!Jinbo Xie debug + do count=1,jall + i = weights_lgr_index_all(count) + + ix = weights_eul_index_all(count,1) + iy = weights_eul_index_all(count,2) + ip = weights_eul_index_all(count,3) + ! + ! convert to 1D indexing of cubed-sphere + ! + ii = (ip-1)*ncube*ncube+(iy-1)*ncube+ix + + wt = weights_all(count,1) + terr_target (i) = terr_target (i) + wt*terr (ii)/area_target(i) +!!Jinbo Xie debug +!#if 0 + landfrac_target (i) = landfrac_target (i) + wt*landfrac (ii)/area_target(i) + landm_coslat_target(i) = landm_coslat_target(i) + wt*landm_coslat(ii)/area_target(i) + sgh30_target (i) = sgh30_target (i) + wt*sgh30 (ii)/area_target(i) +!#endif + tmp = tmp+wt*terr(ii) + end do +!!Jinbo Xie debug +!#endif + + + write(*,*) "tmp", tmp + WRITE(*,*) "max difference between target grid area and remapping software area",& + MAXVAL(target_area-area_target) + + do count=1,ntarget + if (terr_target(count)>8848.0) then + ! + ! max height is higher than Mount Everest + ! + write(*,*) "FATAL error: max height is higher than Mount Everest!" + write(*,*) "terr_target",count,terr_target(count) + write(*,*) "(lon,lat) locations of vertices of cell with excessive max height::" + do i=1,ncorner + write(*,*) target_corner_lon(i,count),target_corner_lat(i,count) + end do + STOP + else if (terr_target(count)<-423.0) then + ! + ! min height is lower than Dead Sea + ! + write(*,*) "FATAL error: min height is lower than Dead Sea!" + write(*,*) "terr_target",count,terr_target(count) + write(*,*) "(lon,lat) locations of vertices of cell with excessive min height::" + do i=1,ncorner + write(*,*) target_corner_lon(i,count),target_corner_lat(i,count) + end do + STOP + else + + end if + end do + WRITE(*,*) "Elevation data passed min/max consistency check!" + WRITE(*,*) + + WRITE(*,*) "min/max of unsmoothed terr_target : ",MINVAL(terr_target ),MAXVAL(terr_target ) + WRITE(*,*) "min/max of landfrac_target : ",MINVAL(landfrac_target),MAXVAL(landfrac_target) + WRITE(*,*) "min/max of landm_coslat_target : ",& + MINVAL(landm_coslat_target),MAXVAL(landm_coslat_target) + WRITE(*,*) "min/max of var30_target : ",MINVAL(sgh30_target ),MAXVAL(sgh30_target ) + ! + ! compute mean height (globally) of topography about sea-level for target grid unfiltered elevation + ! + vol_target_un = 0.0 + area_target_total = 0.0 + DO i=1,ntarget + area_target_total = area_target_total+area_target(i) + vol_target_un = vol_target_un+terr_target(i)*area_target(i) + END DO + WRITE(*,*) "mean height (globally) of topography about sea-level for target grid unfiltered elevation",& + vol_target_un/area_target_total + + ! + ! diagnostics + ! + vol_source = 0.0 + allocate ( dA(ncube,ncube),stat=alloc_error ) + CALL EquiangularAllAreas(ncube, dA) + DO jp=1,6 + DO jy=1,ncube + DO jx=1,ncube + ii = (jp-1)*ncube*ncube+(jy-1)*ncube+jx + vol_source = vol_source+terr(ii)*dA(jx,jy) + END DO + END DO + END DO + WRITE(*,*) "volume of input cubed-sphere terrain :",vol_source + WRITE(*,*) "average elevation of input cubed-sphere terrain:",vol_source/(4.0*pi) + + DEALLOCATE(dA) + ! + ! + ! + allocate (sgh_target(ntarget),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for sgh_target' + stop + end if + ! + ! compute variance with respect to cubed-sphere data + ! + WRITE(*,*) "compute variance with respect to 3km cubed-sphere data: SGH" + + IF (lsmooth_terr) THEN + WRITE(*,*) "smoothing PHIS" + IF (lexternal_smooth_terr) THEN + WRITE(*,*) "using externally generated smoothed topography" + + status = nf_open(trim(smoothed_topography_file), 0, ncid) + IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR(STATUS) + status = nf_close(ncid) + !status = nf_open('phis-smooth.nc', 0, ncid) + !IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR(STATUS) + ! + IF (.NOT.ltarget_latlon) THEN + ! + !********************************************************* + ! + ! read in smoothed topography + ! + !********************************************************* + ! + status = NF_INQ_DIMID (ncid, 'ncol', ntarget_id ) + status = NF_INQ_DIMLEN(ncid, ntarget_id , ntarget_smooth) + IF (ntarget.NE.ntarget_smooth) THEN + WRITE(*,*) "mismatch in smoothed data-set and target grid specification" + WRITE(*,*) ntarget, ntarget_smooth + STOP + END IF + status = NF_INQ_VARID(ncid, 'PHIS', phisid) + ! + ! overwrite terr_target with smoothed version + ! + status = NF_GET_VAR_DOUBLE(ncid, phisid,terr_target) + terr_target = terr_target/9.80616 + ELSE + ! + ! read in smoothed lat-lon topography + ! + status = NF_INQ_DIMID(ncid, 'lon', ntarget_id) + status = NF_INQ_DIMLEN(ncid, ntarget_id, nlon_smooth) + status = NF_INQ_DIMID(ncid, 'lat', ntarget_id) + status = NF_INQ_DIMLEN(ncid, ntarget_id, nlat_smooth) + IF (nlon.NE.nlon_smooth.OR.nlat.NE.nlat_smooth) THEN + WRITE(*,*) "smoothed topography dimensions do not match target grid dimensions" + WRITE(*,*) "target grid : nlon ,nlat =",nlon,nlat + WRITE(*,*) "smoothed topo: nlon_smooth,nlat_smooth =",nlon_smooth,nlat_smooth + STOP + END IF + ALLOCATE(terr_smooth(nlon_smooth,nlat_smooth),stat=alloc_error) + status = NF_INQ_VARID(ncid, 'PHIS', phisid) + status = NF_GET_VAR_DOUBLE(ncid, phisid,terr_smooth) + ! + ! overwrite terr_target with smoothed version + ! + ii=1 + DO j=1,nlat + DO i=1,nlon + terr_target(ii) = terr_smooth(i,j)/9.80616 + ii=ii+1 + END DO + END DO + DEALLOCATE(terr_smooth) + END IF + ELSE + WRITE(*,*) "unstested software - uncomment this line of you know what you are doing!" + STOP + ! + !***************************************************** + ! + ! smoothing topography internally + ! + !***************************************************** + ! + WRITE(*,*) "internally smoothing orography" + ! CALL smooth(terr_target,ntarget,target_corner_lon,target_corner_lat) + ! + ! smooth topography internally + ! + ncoarse = n/(factor*factor) + ! + ! + ! + ncube_coarse = ncube/factor + WRITE(*,*) "resolution of coarse grid", 90.0/ncube_coarse + allocate ( terr_coarse(ncoarse),stat=alloc_error ) + if( alloc_error /= 0 ) then + print*,'Program could not allocate space for landfrac' + stop + end if + WRITE(*,*) "coarsening" + allocate ( dA_coarse(ncube_coarse,ncube_coarse),stat=alloc_error ) + CALL coarsen(terr,terr_coarse,factor,n,dA_coarse) + ! + ! + ! + vol_tmp = 0.0 + DO jp=1,6 + DO jy=1,ncube_coarse + DO jx=1,ncube_coarse + ii = (jp-1)*ncube_coarse*ncube_coarse+(jy-1)*ncube_coarse+jx + vol_tmp = vol_tmp+terr_coarse(ii)*dA_coarse(jx,jy) + END DO + END DO + END DO + WRITE(*,*) "volume of coarsened cubed-sphere terrain :",vol_source + WRITE(*,*) "difference between coarsened cubed-sphere data and input cubed-sphere data",& + vol_tmp-vol_source + + + + WRITE(*,*) "done coarsening" + + nreconstruction = 1 + IF (norder>1) THEN + IF (norder == 2) THEN + nreconstruction = 3 + ELSEIF (norder == 3) THEN + nreconstruction = 6 + END IF + ALLOCATE(recons (nreconstruction, ncoarse), STAT=status) + ALLOCATE(centroids(nreconstruction, ncoarse), STAT=status) + CALL get_reconstruction(terr_coarse,norder, nmono, recons, npd,da_coarse,& + ncube_coarse+1,nreconstruction,centroids) + SELECT CASE (nmono) + CASE (0) + WRITE(*,*) "coarse grid reconstructions are not filtered with shape-preesrving filter" + CASE (1) + WRITE(*,*) "coarse grid reconstructions are filtered with shape-preserving filter" + CASE DEFAULT + WRITE(*,*) "nmono out of range: ",nmono + STOP + END SELECT + SELECT CASE (0) + CASE (0) + WRITE(*,*) "coarse grid reconstructions are not filtered with positive definite filter" + CASE (1) + WRITE(*,*) "coarse grid reconstructions filtered with positive definite filter" + CASE DEFAULT + WRITE(*,*) "npd out of range: ",npd + STOP + END SELECT + END IF + + jall_coarse = (ncube*ncube*12) !anticipated number of weights + jmax_segments_coarse = jmax_segments!/factor ! + WRITE(*,*) "anticipated",jall_coarse + allocate (weights_all_coarse(jall_coarse,nreconstruction),stat=alloc_error ) + allocate (weights_eul_index_all_coarse(jall_coarse,3),stat=alloc_error ) + allocate (weights_lgr_index_all_coarse(jall_coarse),stat=alloc_error ) + ! + ! + ! + CALL overlap_weights(weights_lgr_index_all_coarse,weights_eul_index_all_coarse,weights_all_coarse,& + jall_coarse,ncube_coarse,ngauss,ntarget,ncorner,jmax_segments_coarse,target_corner_lon,& + target_corner_lat,nreconstruction) + + WRITE(*,*) "MIN/MAX of area-weight [0:1]: ",& + MINVAL(weights_all_coarse(:,1)),MAXVAL(weights_all_coarse(:,1)) + ! + ! compute new weights + ! + + ! + ! do mapping + ! + terr_target = 0.0 + tmp = 0.0 + allocate ( area_target_coarse(ntarget),stat=alloc_error) + all_weights = 0.0 + area_target_coarse = 0.0 + do count=1,jall_coarse + i = weights_lgr_index_all_coarse(count) + wt = weights_all_coarse(count,1) + area_target_coarse (i) = area_target_coarse(i) + wt + all_weights = all_weights+wt + end do + WRITE(*,*) "sum of all weights (coarse to target) minus area of sphere : ",all_weights-4.0*pi + WRITE(*,*) "MIN/MAX of area_target_coarse [0:1]:",& + MINVAL(area_target_coarse),MAXVAL(area_target_coarse) + IF (norder==1) THEN + do count=1,jall_coarse + i = weights_lgr_index_all_coarse(count) + + ix = weights_eul_index_all_coarse(count,1) + iy = weights_eul_index_all_coarse(count,2) + ip = weights_eul_index_all_coarse(count,3) + ! + ! convert to 1D indexing of cubed-sphere + ! + ii = (ip-1)*ncube_coarse*ncube_coarse+(iy-1)*ncube_coarse+ix + + wt = weights_all_coarse(count,1) + + terr_target(i) = terr_target(i) + wt*terr_coarse(ii)/area_target_coarse(i) + tmp = tmp+wt*terr_coarse(ii) + end do + ELSE IF (norder==2) THEN + do count=1,jall_coarse + i = weights_lgr_index_all_coarse(count) + IF (i>jall_coarse.OR.i<1) THEN + WRITE(*,*) i,jall_coarse + STOP + END IF + ix = weights_eul_index_all_coarse(count,1) + iy = weights_eul_index_all_coarse(count,2) + ip = weights_eul_index_all_coarse(count,3) + ! + ! convert to 1D indexing of cubed-sphere + ! + ii = (ip-1)*ncube_coarse*ncube_coarse+(iy-1)*ncube_coarse+ix + + terr_target(i) = terr_target(i) + (weights_all_coarse(count,1)*(& + ! + ! all constant terms + ! + terr_coarse(ii) & + - recons(1,ii)*centroids(1,ii) & + - recons(2,ii)*centroids(2,ii) & + ! + ! + recons(3,ii)*(2.0*centroids(1,ii)**2-centroids(3,ii))& + ! + recons(4,ii)*(2.0*centroids(2,ii)**2-centroids(4,ii))& + ! + ! + recons(5,ii)*(2.0*centroids(1,ii)*centroids(2,ii)-centroids(5,ii))& + )+& + ! + ! linear terms + ! + weights_all_coarse(count,2)*(& + + recons(1,ii)& + + ! - recons(3,ii)*2.0*centroids(1,ii)& + ! - recons(5,ii)* centroids(2,ii)& + )+& + ! + weights_all_coarse(count,3)*(& + recons(2,ii)& + ! + ! - recons(4,ii)*2.0*centroids(2,ii)& + ! - recons(5,ii)* centroids(1,ii)& + )& + ! + ! quadratic terms + ! + ! weights_all_coarse(count,4)*recons(3,ii)+& + ! weights_all_coarse(count,5)*recons(4,ii)+& + ! weights_all_coarse(count,6)*recons(5,ii) + )/area_target_coarse(i) + end do + DEALLOCATE(centroids) + DEALLOCATE(recons) + DEALLOCATE(weights_all_coarse) + + ELSE IF (norder==3) THEN + ! recons(4,:) = 0.0 + ! recons(5,:) = 0.0 + do count=1,jall_coarse + i = weights_lgr_index_all_coarse(count) + IF (i>jall_coarse.OR.i<1) THEN + WRITE(*,*) i,jall_coarse + STOP + END IF + ix = weights_eul_index_all_coarse(count,1) + iy = weights_eul_index_all_coarse(count,2) + ip = weights_eul_index_all_coarse(count,3) + ! + ! convert to 1D indexing of cubed-sphere + ! + ii = (ip-1)*ncube_coarse*ncube_coarse+(iy-1)*ncube_coarse+ix + + ! terr_target(i) = terr_target(i) + wt*terr_coarse(ii)/area_target_coarse(i) + + ! WRITE(*,*) count,area_target_coarse(i) + ! terr_target(i) = terr_target(i) + area_target_coarse(i) + ! + terr_target(i) = terr_target(i) + (weights_all_coarse(count,1)*(& + + + ! centroids(5,ii))/area_target_coarse(i)) + ! centroids(1,ii)/area_target_coarse(i)) + ! /area_target_coarse(i)) + + + + + ! + ! all constant terms + ! + terr_coarse(ii) & + - recons(1,ii)*centroids(1,ii) & + - recons(2,ii)*centroids(2,ii) & + ! + + recons(3,ii)*(2.0*centroids(1,ii)**2-centroids(3,ii))& + + recons(4,ii)*(2.0*centroids(2,ii)**2-centroids(4,ii))& + ! + + recons(5,ii)*(2.0*centroids(1,ii)*centroids(2,ii)-centroids(5,ii))& + )+& + ! + ! linear terms + ! + weights_all_coarse(count,2)*(& + + recons(1,ii)& + + - recons(3,ii)*2.0*centroids(1,ii)& + - recons(5,ii)* centroids(2,ii)& + )+& + ! + weights_all_coarse(count,3)*(& + recons(2,ii)& + ! + - recons(4,ii)*2.0*centroids(2,ii)& + - recons(5,ii)* centroids(1,ii)& + )+& + ! + ! quadratic terms + ! + weights_all_coarse(count,4)*recons(3,ii)+& + weights_all_coarse(count,5)*recons(4,ii)+& + weights_all_coarse(count,6)*recons(5,ii))/area_target_coarse(i) + end do + DEALLOCATE(centroids) + DEALLOCATE(recons) + DEALLOCATE(weights_all_coarse) + END IF + DEALLOCATE(area_target_coarse) + WRITE(*,*) "done smoothing" + END IF + ! + ! compute mean height (globally) of topography about sea-level for target grid filtered elevation + ! + vol_target = 0.0 + DO i=1,ntarget + vol_target = vol_target+terr_target(i)*area_target(i) + ! if (ABS(area_target(i)-area_target_coarse(i))>0.000001) THEN + ! WRITE(*,*) "xxx",area_target(i),area_target_coarse(i),area_target(i)-area_target_coarse(i) + ! STOP + ! END IF + END DO + WRITE(*,*) "mean height (globally) of topography about sea-level for target grid filtered elevation",& + vol_target/area_target_total + WRITE(*,*) "percentage change in mean height between filtered and unfiltered elevations",& + 100.0*(vol_target-vol_target_un)/vol_target_un + WRITE(*,*) "percentage change in mean height between input cubed-sphere and unfiltered elevations",& + 100.0*(vol_source-vol_target_un)/vol_source + + END IF + ! + ! Done internal smoothing + ! + WRITE(*,*) "min/max of terr_target : ",MINVAL(terr_target),MAXVAL(terr_target) + + if (lzero_out_ocean_point_phis) then + WRITE(*,*) "if ocean mask PHIS=0.0" + end if + + + sgh_target=0.0 + do count=1,jall + i = weights_lgr_index_all(count)!! + ! + ix = weights_eul_index_all(count,1) + iy = weights_eul_index_all(count,2) + ip = weights_eul_index_all(count,3) + ! + ! convert to 1D indexing of cubed-sphere + ! + ii = (ip-1)*ncube*ncube+(iy-1)*ncube+ix! + + wt = weights_all(count,1) + + if (lzero_out_ocean_point_phis.AND.landfrac_target(i).lt.0.01_r8) then + terr_target(i) = 0.0_r8 !5*terr_target(i) + end if + sgh_target(i) = sgh_target(i)+wt*((terr_target(i)-terr(ii))**2)/area_target(i) + end do + + + + + ! + ! zero out small values + ! + DO i=1,ntarget + IF (landfrac_target(i)<.001_r8) landfrac_target(i) = 0.0 + IF (sgh_target(i)<0.5) sgh_target(i) = 0.0 + IF (sgh30_target(i)<0.5) sgh30_target(i) = 0.0 + END DO + sgh_target = SQRT(sgh_target) + sgh30_target = SQRT(sgh30_target) + +!============Jinbo Xie========== +!for centroid of mass +!wt is useful proxy for dA +!#if 0 +!!Jinbo Xie debug +print*,"cal oa" +allocate(oa_target(ntarget,nvar_dirOA),stat=alloc_error) +call OAdir(terr,ntarget,ncube,n,nvar_dirOA,jall,weights_lgr_index_all,weights_eul_index_all(:,1),weights_eul_index_all(:,2),weights_eul_index_all(:,3),weights_all,landfrac_target,target_center_lon,target_center_lat,lon_terr,lat_terr,area_target,oa_target)!OAx,OAy) +!call OAorig(terr,ntarget,ncube,n,jall,weights_lgr_index_all,weights_eul_index_all(:,1),weights_eul_index_all(:,2),weights_eul_index_all(:,3),weights_all,landfrac_target,lon_terr,lat_terr,area_target,oa_target) +!#endif +!============Jinbo Xie par========== +!Jinbo Xie +!OC + print*,"cal oc" + allocate(oc_target(ntarget),stat=alloc_error) + oc_target=0.0_r8 + call OC(terr,ntarget,ncube,n,jall,weights_lgr_index_all,weights_eul_index_all(:,1),weights_eul_index_all(:,2),weights_eul_index_all(:,3),weights_all,landfrac_target,area_target,sgh_target,terr_target,oc_target) + +!stop +!#if 0 +!#endif +!OL + print*,"cal ol" + allocate(ol_target(ntarget,nvar_dirOL),stat=alloc_error) + ol_target=0.0_r8 + !call OLorig(terr,ntarget,ncube,n,jall,weights_lgr_index_all,weights_eul_index_all(:,1),weights_eul_index_all(:,2),weights_eul_index_all(:,3),weights_all,landfrac_target,lon_terr,lat_terr,area_target,sgh_target,target_center_lat,target_center_lon,target_corner_lat_deg,target_corner_lon_deg,ol_target) +!!Jinbo Xie debug +!#endif + allocate(indexb(ntarget),stat=alloc_error) + indexb=0.0_r8 + do count=1,jall + i = weights_lgr_index_all(count) + indexb(i)=indexb(i)+1 + enddo + allocate(terrout(4,ntarget,maxval(indexb)),stat=alloc_error) + allocate(dxy(ntarget,nvar_dirOL),stat=alloc_error) + !call OLdir(terr,ntarget,ncube,n,jall,nlon,nlat,maxval(indexb),nvar_dirOL,weights_lgr_index_all,weights_eul_index_all(:,1),weights_eul_index_all(:,2),weights_eul_index_all(:,3),weights_all,landfrac_target,target_center_lon,target_center_lat,lon_terr,lat_terr,sgh_target,ol_target,terrout,dxy) + !call OLdir(terr,ntarget,ncube,n,jall,nlon,nlat,maxval(indexb),nvar_dirOL,weights_lgr_index_all,weights_eul_index_all(:,1),weights_eul_index_all(:,2),weights_eul_index_all(:,3),weights_all,landfrac_target,target_center_lon,target_center_lat,lon_terr,lat_terr,sgh_target,ol_target,terrout) + !call OLorig(terr,ntarget,ncube,n,jall,weights_lgr_index_all,weights_eul_index_all(:,1),weights_eul_index_all(:,2),weights_eul_index_all(:,3),weights_all,landfrac_target,lon_terr,lat_terr,area_target,sgh_target,target_center_lat,target_center_lon,target_corner_lat_deg,target_corner_lon_deg,ol_target) + !call OLdir(terr,ntarget,ncube,n,jall,nlon,nlat,maxval(indexb),nvar_dirOL,weights_lgr_index_all,weights_eul_index_all(:,1),weights_eul_index_all(:,2),weights_eul_index_all(:,3),weights_all,landfrac_target,target_center_lon,target_center_lat,lon_terr,lat_terr,sgh_target,area_target,ol_target,terrout,dxy) + call OLdir(terr,ntarget,ncube,n,jall,nlon,nlat,maxval(indexb),nvar_dirOL,weights_lgr_index_all,weights_eul_index_all(:,1),weights_eul_index_all(:,2),weights_eul_index_all(:,3),weights_all,landfrac_target,target_center_lon,target_center_lat,target_corner_lon_deg,target_corner_lat_deg,lon_terr,lat_terr,sgh_target,area_target,ol_target,terrout,dxy) + !do i=1,10!180 + !print*,"OLdir Jinbo Xie",minval(ol_target(:,i)),maxval(ol_target(:,i)) + !enddo + !stop +!!Jinbo Xie debug +!#endif +!#endif +!========Jinbo Xie par========= + + + WRITE(*,*) "min/max of sgh_target : ",MINVAL(sgh_target),MAXVAL(sgh_target) + WRITE(*,*) "min/max of sgh30_target : ",MINVAL(sgh30_target),MAXVAL(sgh30_target) + + DEALLOCATE(terr,weights_all,weights_eul_index_all,landfrac,landm_coslat) + + + + IF (ltarget_latlon) THEN +#if 0 + CALL wrtncdf_rll(nlon,nlat,lpole,ntarget,terr_target,landfrac_target,sgh_target,sgh30_target,& + landm_coslat_target,target_center_lon,target_center_lat,.true.) +#endif +!========Jinbo Xie========== +print*,"output rll" + CALL wrtncdf_rll(nlon,nlat,nvar_dirOA,nvar_dirOL,maxval(indexb),lpole,ntarget,terr_target,landfrac_target,sgh_target,sgh30_target, oc_target,oa_target,ol_target,terrout,dxy,& + landm_coslat_target,target_center_lon,target_center_lat,.false.,output_topography_file) +!========Jinbo Xie========== + + ELSE +#if 0 + CALL wrtncdf_unstructured(ntarget,terr_target,landfrac_target,sgh_target,sgh30_target,& + landm_coslat_target,target_center_lon,target_center_lat) +#endif +!========Jinbo Xie========== + print*,"output unstructure" + CALL wrtncdf_unstructured(nvar_dirOA,nvar_dirOL,maxval(indexb),ntarget,terr_target,landfrac_target,sgh_target,sgh30_target,oc_target,oa_target,ol_target,terrout,dxy,landm_coslat_target,target_center_lon,target_center_lat,output_topography_file) +!========Jinbo Xie========== + END IF + + DEALLOCATE(terr_target,landfrac_target,sgh30_target,sgh_target,landm_coslat_target) + +!====Jinbo Xie==== +DEALLOCATE(oc_target) +!====Jinbo Xie==== + +end program convterr + +! +! +! +#if 0 +subroutine wrtncdf_unstructured(n,terr,landfrac,sgh,sgh30,landm_coslat,lon,lat) +#endif +subroutine wrtncdf_unstructured(nvar_dirOA,nvar_dirOL,indexb,n,terr,landfrac,sgh,sgh30,oc_in,oa_in,ol_in,terrout,dxy_in,landm_coslat,lon,lat,output) + + use shr_kind_mod, only: r8 => shr_kind_r8 + implicit none + +# include + + ! + ! Dummy arguments + ! + integer, intent(in) :: n + real(r8),dimension(n) , intent(in) :: terr, landfrac,sgh,sgh30,lon, lat, landm_coslat + ! + ! Local variables + ! + character (len=512) :: fout ! NetCDF output file + + integer :: foutid ! Output file id + integer :: lonid, lonvid + integer :: latid, latvid + integer :: terrid,nid + integer :: terrdim,landfracid,sghid,sgh30id,landm_coslatid + integer :: status ! return value for error control of netcdf routin + integer :: i,j + integer, dimension(2) :: nc_lat_vid,nc_lon_vid + character (len=8) :: datestring + integer :: nc_gridcorn_id, lat_vid, lon_vid + + real(r8), parameter :: fillvalue = 1.d36 + !=====Jinbo Xie======== + integer, intent(in) :: nvar_dirOA,nvar_dirOL,indexb + character(len=512) :: output + !Jinbo Xie add direction + !=====Jinbo Xie======== + integer :: ocid,varid,var2id,indexbid,terroutid(4) + integer :: oaid,olid,dxyid + integer :: oa1id,oa2id,oa3id,oa4id + integer :: ol1id,ol2id,ol3id,ol4id + !======Jinbo Xie======= + integer, dimension(2) :: ocdim + integer, dimension(3) :: oadim,oldim,terroutdim + !======Jinbo Xie========= + real(r8),dimension(n) , intent(in) :: oc_in + real(r8),dimension(n,nvar_dirOA) , intent(in) :: oa_in + real(r8),dimension(n,nvar_dirOL) , intent(in) :: ol_in + real(r8),dimension(4,n,indexb),intent(in) :: terrout + real(r8),dimension(n,nvar_dirOL),intent(in) :: dxy_in + character*20,dimension(4) :: terroutchar + !!=======Jinbo Xie========= + real(r8),dimension(n) :: oc + real(r8),dimension(n,nvar_dirOA) :: oa + real(r8),dimension(n,nvar_dirOL) :: ol + real(r8),dimension(n,nvar_dirOL) :: dxy + character*20 :: numb + !!======Jinbo Xie======= + write(numb,"(i0.1)") nvar_dirOL + print*,"dir number", nvar_dirOL + !fout='final-'//adjustl(trim(numb))//'.nc' + fout=output + !!======Jinbo Xie======== + !print*,"Jinbo Xie shape(oc_in),shape(oc)",shape(oc_in),shape(oc) + oc=oc_in + oa=oa_in + ol=ol_in + dxy=dxy_in + !Jinbo Xie debug + !!======Jinbo Xie======== + ! + ! Create NetCDF file for output + ! + print *,"Create NetCDF file for output" + status = nf_create (fout, NF_64BIT_OFFSET , foutid) + if (status .ne. NF_NOERR) call handle_err(status) + ! + ! Create dimensions for output + ! + status = nf_def_dim (foutid, 'ncol', n, nid) + if (status .ne. NF_NOERR) call handle_err(status) + + !!====Jinbo Xie======== + status = nf_def_dim (foutid, 'nvar_dirOA', nvar_dirOA, varid) + if (status .ne. NF_NOERR) call handle_err(status) + status = nf_def_dim (foutid, 'nvar_dirOL', nvar_dirOL, var2id) + if (status .ne. NF_NOERR) call handle_err(status) + + !Jinbo Xie debug + !status = nf_def_dim (foutid, 'indexb',23, indexbid) + status = nf_def_dim (foutid, 'indexb', indexb, indexbid) + !Jinbo Xie debug + if (status .ne. NF_NOERR) call handle_err(status) + !!=====Jinbo Xie===== + ! + ! Create variable for output + ! + print *,"Create variable for output" + status = nf_def_var (foutid,'PHIS', NF_DOUBLE, 1, nid, terrid) + if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_def_var (foutid,'LANDFRAC', NF_DOUBLE, 1, nid, landfracid) + if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_def_var (foutid,'SGH', NF_DOUBLE, 1, nid, sghid) + if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_def_var (foutid,'SGH30', NF_DOUBLE, 1, nid, sgh30id) + if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_def_var (foutid,'LANDM_COSLAT', NF_DOUBLE, 1, nid, landm_coslatid) + if (status .ne. NF_NOERR) call handle_err(status) + ! + status = nf_def_var (foutid,'lat', NF_DOUBLE, 1, nid, latvid) + if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_def_var (foutid,'lon', NF_DOUBLE, 1, nid, lonvid) + if (status .ne. NF_NOERR) call handle_err(status) + + !!========Jinbo Xie======== + status = nf_def_var (foutid,'OC', NF_DOUBLE, 1, nid, ocid) + oadim(1)=nid + oadim(2)=varid + status = nf_def_var (foutid,'OA', NF_DOUBLE, 2, oadim, oaid) + oldim(1)=nid + oldim(2)=var2id + status = nf_def_var (foutid,'OL', NF_DOUBLE, 2, oldim, olid) +#if 0 + terroutdim(1)=nid + terroutdim(2)=indexbid + !name + terroutchar(1)="terr" + terroutchar(2)="terrx" + terroutchar(3)="terry" + terroutchar(4)="wt" + do i=1,4 + status = nf_def_var (foutid, terroutchar(i), NF_DOUBLE, 2, & + terroutdim, terroutid(i)) + enddo + !dxy + status = nf_def_var (foutid,'dxy', NF_DOUBLE, 2, oldim, dxyid) +#endif + !!========Jinbo Xie========== + ! + ! Create attributes for output variables + ! + status = nf_put_att_text (foutid,terrid,'long_name', 21, 'surface geopotential') + status = nf_put_att_text (foutid,terrid,'units', 5, 'm2/s2') + status = nf_put_att_double (foutid, terrid, 'missing_value', nf_double, 1, fillvalue) + status = nf_put_att_double (foutid, terrid, '_FillValue' , nf_double, 1, fillvalue) + ! status = nf_put_att_text (foutid,terrid,'filter', 35, 'area averaged from USGS 30-sec data') + + status = nf_put_att_double (foutid, sghid, 'missing_value', nf_double, 1, fillvalue) + status = nf_put_att_double (foutid, sghid, '_FillValue' , nf_double, 1, fillvalue) + status = nf_put_att_text (foutid, sghid, 'long_name' , 48, & + 'standard deviation of 3km cubed-sphere elevation and target grid elevation') + status = nf_put_att_text (foutid, sghid, 'units' , 1, 'm') + ! status = nf_put_att_text (foutid, sghid, 'filter' , 4, 'none') + + status = nf_put_att_double (foutid, sgh30id, 'missing_value', nf_double, 1, fillvalue) + status = nf_put_att_double (foutid, sgh30id, '_FillValue' , nf_double, 1, fillvalue) + status = nf_put_att_text (foutid, sgh30id, 'long_name' , 49, & + 'standard deviation of 30s elevation from 3km cubed-sphere cell average height') + status = nf_put_att_text (foutid, sgh30id, 'units' , 1, 'm') + ! status = nf_put_att_text (foutid, sgh30id, 'filter' , 4, 'none') + + status = nf_put_att_double (foutid, landm_coslatid, 'missing_value', nf_double, 1, fillvalue) + status = nf_put_att_double (foutid, landm_coslatid, '_FillValue' , nf_double, 1, fillvalue) + status = nf_put_att_text (foutid, landm_coslatid, 'long_name' , 23, 'smoothed land fraction') + status = nf_put_att_text (foutid, landm_coslatid, 'filter' , 4, 'none') + + status = nf_put_att_double (foutid, landfracid, 'missing_value', nf_double, 1, fillvalue) + status = nf_put_att_double (foutid, landfracid, '_FillValue' , nf_double, 1, fillvalue) + status = nf_put_att_text (foutid, landfracid, 'long_name', 21, 'gridbox land fraction') + ! status = nf_put_att_text (foutid, landfracid, 'filter', 40, 'area averaged from 30-sec USGS raw data') + + + status = nf_put_att_text (foutid,latvid,'long_name', 8, 'latitude') + if (status .ne. NF_NOERR) call handle_err(status) + status = nf_put_att_text (foutid,latvid,'units', 13, 'degrees_north') + if (status .ne. NF_NOERR) call handle_err(status) + ! status = nf_put_att_text (foutid,latvid,'units', 21, 'cell center locations') + ! if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_put_att_text (foutid,lonvid,'long_name', 9, 'longitude') + if (status .ne. NF_NOERR) call handle_err(status) + status = nf_put_att_text (foutid,lonvid,'units', 12, 'degrees_east') + if (status .ne. NF_NOERR) call handle_err(status) + ! status = nf_put_att_text (foutid,lonvid,'units' , 21, 'cell center locations') + ! if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_put_att_text (foutid,NF_GLOBAL,'source', 50, 'USGS 30-sec dataset binned to ncube3000 (cube-sphere) grid') + if (status .ne. NF_NOERR) call handle_err(status) + status = nf_put_att_text (foutid,NF_GLOBAL,'title', 24, '30-second USGS topo data') + if (status .ne. NF_NOERR) call handle_err(status) + call DATE_AND_TIME(DATE=datestring) + status = nf_put_att_text (foutid,NF_GLOBAL,'history',25, 'Written on date: ' // datestring ) + if (status .ne. NF_NOERR) call handle_err(status) + !!======Jinbo Xie============ + status = nf_put_att_text (foutid,oaid,'note', 40, '(2)+1 in nvar_dirOA to avoid bug in io') +#if 0 + do i=1,4 + status = nf_put_att_double (foutid, terroutid(i),& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, terroutid(i),& + '_FillValue' , nf_double, 1,fillvalue) + enddo +#endif + !!======Jinbo Xie============ + + ! + ! End define mode for output file + ! + status = nf_enddef (foutid) + if (status .ne. NF_NOERR) call handle_err(status) + ! + ! Write variable for output + ! + + !!==========Jinbo Xie============ + print*,"writing oc data",MINVAL(oc),MAXVAL(oc) + status = nf_put_var_double (foutid, ocid, oc) + if (status .ne. NF_NOERR) call handle_err(status) + !oa,ol + print*,"writing oa data",MINVAL(oa),MAXVAL(oa) + status = nf_put_var_double (foutid, oaid, oa) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"writing ol data",MINVAL(ol),MAXVAL(ol) + status = nf_put_var_double (foutid, olid, ol) + !========Jinbo Xie======== + !=========== + if (status .ne. NF_NOERR) call handle_err(status) +#if 0 + do i=1,4 + status = nf_put_att_double (foutid, terroutid(i),& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, terroutid(i),& + '_FillValue' , nf_double, 1,fillvalue) + print*,"writing"//terroutchar(i)//" data",& + MINVAL(terrout(i,:,:)),MAXVAL(terrout(i,:,:)) + status = nf_put_var_double (foutid, terroutid(i), terrout(i,:,:)) + if (status .ne. NF_NOERR) call handle_err(status) + enddo +!#endif +!#if 0 + print*,"writing dxy data",MINVAL(dxy),MAXVAL(dxy) + status = nf_put_var_double (foutid, dxyid, dxy) + if (status .ne. NF_NOERR) call handle_err(status) +#endif +!========Jinbo Xie======== + + print*,"writing terrain data",MINVAL(terr),MAXVAL(terr) + status = nf_put_var_double (foutid, terrid, terr*9.80616) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing terrain data" + + print*,"writing landfrac data",MINVAL(landfrac),MAXVAL(landfrac) + status = nf_put_var_double (foutid, landfracid, landfrac) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing landfrac data" + + print*,"writing sgh data",MINVAL(sgh),MAXVAL(sgh) + status = nf_put_var_double (foutid, sghid, sgh) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing sgh data" + + print*,"writing sgh30 data",MINVAL(sgh30),MAXVAL(sgh30) + status = nf_put_var_double (foutid, sgh30id, sgh30) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing sgh30 data" + + print*,"writing landm_coslat data",MINVAL(landm_coslat),MAXVAL(landm_coslat) + status = nf_put_var_double (foutid, landm_coslatid, landm_coslat) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing sgh30 data" + ! + print*,"writing lat data" + status = nf_put_var_double (foutid, latvid, lat) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing lat data" + + print*,"writing lon data" + status = nf_put_var_double (foutid, lonvid, lon) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing lon data" + ! + ! Close output file + ! + print *,"close file" + status = nf_close (foutid) + if (status .ne. NF_NOERR) call handle_err(status) +end subroutine wrtncdf_unstructured +! +!************************************************************** +! +! if target grid is lat-lon output structured +! +!************************************************************** +! + +!=======Jinbo Xie========= +#if 0 +subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,landm_coslat_in,lon,lat,lprepare_fv_smoothing_routine) +#endif +!=======Jinbo Xie========= +subroutine wrtncdf_rll(nlon,nlat,nvar_dirOA,nvar_dirOL,indexb,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,oc_in,oa_in,ol_in,terrout,dxy_in,landm_coslat_in,lon,lat,lprepare_fv_smoothing_routine,output) +!=======Jinbo Xie========= + + use shr_kind_mod, only: r8 => shr_kind_r8 + implicit none + +# include + + ! + ! Dummy arguments + ! + integer, intent(in) :: n,nlon,nlat,nvar_dirOA,nvar_dirOL,indexb + !Jinbo Xie add direction + ! + ! lprepare_fv_smoothing_routine is to make a NetCDF file that can be used with the CAM-FV smoothing software + ! + logical , intent(in) :: lpole,lprepare_fv_smoothing_routine + real(r8),dimension(n) , intent(in) :: terr_in, landfrac_in,sgh_in,sgh30_in,lon, lat, landm_coslat_in + +!=======Jinbo Xie========= + real(r8),dimension(n) , intent(in) :: oc_in + real(r8),dimension(n,nvar_dirOA) , intent(in) :: oa_in + real(r8),dimension(n,nvar_dirOL) , intent(in) :: ol_in + real(r8),dimension(4,n,indexb),intent(in) :: terrout + real(r8),dimension(n,nvar_dirOL),intent(in) :: dxy_in + character*20,dimension(4) :: terroutchar + character(len=512),intent(in) :: output +!=======Jinbo Xie========= + + ! + ! Local variables + ! + character (len=512):: fout ! NetCDF output file + integer :: foutid ! Output file id + integer :: lonid, lonvid + integer :: latid, latvid + integer :: terrid,nid +!=====Jinbo Xie======== + integer :: ocid,varid,var2id,indexbid,terroutid(4) + integer :: oaid,olid,dxyid +integer :: oa1id,oa2id,oa3id,oa4id +integer :: ol1id,ol2id,ol3id,ol4id +!=====Jinbo Xie======== + integer :: terrdim,landfracid,sghid,sgh30id,landm_coslatid + integer :: status ! return value for error control of netcdf routin + integer :: i,j + integer, dimension(2) :: nc_lat_vid,nc_lon_vid + character (len=8) :: datestring + integer :: nc_gridcorn_id, lat_vid, lon_vid + real(r8), parameter :: fillvalue = 1.d36 + real(r8) :: ave + + real(r8),dimension(nlon) :: lonar ! longitude array + real(r8),dimension(nlat) :: latar ! latitude array + + integer, dimension(2) :: htopodim,landfdim,sghdim,sgh30dim,landmcoslatdim +!======Jinbo Xie======= +integer, dimension(2) :: ocdim +integer, dimension(3) :: oadim,oldim,terroutdim +!======Jinbo Xie======= + real(r8),dimension(n) :: terr, landfrac,sgh,sgh30,landm_coslat +!======Jinbo Xie======= + real(r8),dimension(n) :: oc + real(r8),dimension(n,nvar_dirOA) :: oa + real(r8),dimension(n,nvar_dirOL) :: ol + real(r8),dimension(n,nvar_dirOL) :: dxy + character*20 :: numb +!======Jinbo Xie======= + +!print*,"nlon nlat n",nlon, nlat, n + IF (nlon*nlat.NE.n) THEN + WRITE(*,*) "inconsistent input for wrtncdf_rll" + STOP + END IF + ! + ! we assume that the unstructured layout of the lat-lon grid is ordered in latitude rows, that is, + ! unstructured index n is given by + ! + ! n = (j-1)*nlon+i + ! + ! where j is latitude index and i longitude index + ! + do i = 1,nlon + lonar(i)= lon(i) + enddo + do j = 1,nlat + latar(j)= lat((j-1)*nlon+1) + enddo + + terr = terr_in + sgh=sgh_in + sgh30 =sgh30_in + landfrac = landfrac_in + landm_coslat = landm_coslat_in + +!====Jinbo Xie====== + oc=oc_in + oa=oa_in + ol=ol_in + dxy=dxy_in +!====Jinbo Xie====== + + + if (lpole) then + write(*,*) "average pole control volume" + ! + ! North pole - terr + ! + ave = 0.0 + do i=1,nlon + ave = ave + terr_in(i) + end do + terr(1:nlon) = ave/DBLE(nlon) + ! + ! South pole + ! + ave = 0.0 + do i=n-(nlon+1),n + ave = ave + terr_in(i) + end do + terr(n-(nlon+1):n) = ave/DBLE(nlon) + + +!=========Jinbo Xie========= + !oc + ! North pole - terr + ave = 0.0 + do i=1,nlon + ave = ave + oc_in(i) + end do + oc(1:nlon) = ave/DBLE(nlon) + ! South pole + ave = 0.0 + do i=n-(nlon+1),n + ave = ave + oc_in(i) + end do + oc(n-(nlon+1):n) = ave/DBLE(nlon) + !oa + ! North pole - terr +do j =1,nvar_dirOA + ave = 0.0 + do i=1,nlon + ave = ave + oa_in(i,j) + end do + oa(1:nlon,j) = ave/DBLE(nlon) + ! South pole + ave = 0.0 + do i=n-(nlon+1),n + ave = ave + oa_in(i,j) + end do + oa(n-(nlon+1):n,j) = ave/DBLE(nlon) +enddo + !ol +!#if 0 +! North pole - terr +do j =1,nvar_dirOL + ave = 0.0 + do i=1,nlon + ave = ave + ol_in(i,j) + end do + ol(1:nlon,j) = ave/DBLE(nlon) + ! South pole + ave = 0.0 + do i=n-(nlon+1),n + ave = ave + ol_in(j,i) + end do + ol(n-(nlon+1):n,j) = ave/DBLE(nlon) +enddo +!#endif +!=========Jinbo Xie========= + + ! + ! North pole - sgh + ! + ave = 0.0 + do i=1,nlon + ave = ave + sgh_in(i) + end do + sgh(1:nlon) = ave/DBLE(nlon) + ! + ! South pole + ! + ave = 0.0 + do i=n-(nlon+1),n + ave = ave + sgh_in(i) + end do + sgh(n-(nlon+1):n) = ave/DBLE(nlon) + + ! + ! North pole - sgh30 + ! + ave = 0.0 + do i=1,nlon + ave = ave + sgh30_in(i) + end do + sgh30(1:nlon) = ave/DBLE(nlon) + ! + ! South pole + ! + ave = 0.0 + do i=n-(nlon+1),n + ave = ave + sgh30_in(i) + end do + sgh30(n-(nlon+1):n) = ave/DBLE(nlon) + + ! + ! North pole - landfrac + ! + ave = 0.0 + do i=1,nlon + ave = ave + landfrac_in(i) + end do + landfrac(1:nlon) = ave/DBLE(nlon) + ! + ! South pole + ! + ave = 0.0 + do i=n-(nlon+1),n + ave = ave + landfrac_in(i) + end do + landfrac(n-(nlon+1):n) = ave/DBLE(nlon) + + ! + ! North pole - landm_coslat + ! + ave = 0.0 + do i=1,nlon + ave = ave + landm_coslat_in(i) + end do + landm_coslat(1:nlon) = ave/DBLE(nlon) + ! + ! South pole + ! + ave = 0.0 + do i=n-(nlon+1),n + ave = ave + landm_coslat_in(i) + end do + landm_coslat(n-(nlon+1):n) = ave/DBLE(nlon) + +!dxy Jinbo Xie + do j=1,4 + ave = 0.0 + do i=n-(nlon+1),n + ave = ave + dxy(j,i) + end do + dxy(j,n-(nlon+1):n) = ave/DBLE(nlon) + enddo +!dxy Jinbo Xie + end if + +!print Jinbo + write(numb,"(i0.1)") nvar_dirOL + print*,"dir number", nvar_dirOL + + + !fout='final-'//adjustl(trim(numb))//'.nc' + fout=output + ! + ! Create NetCDF file for output + ! + print *,"Create NetCDF file for output" + status = nf_create (fout, NF_64BIT_OFFSET , foutid) + if (status .ne. NF_NOERR) call handle_err(status) + ! + ! Create dimensions for output + ! + print *,"Create dimensions for output" + status = nf_def_dim (foutid, 'lon', nlon, lonid) + if (status .ne. NF_NOERR) call handle_err(status) + status = nf_def_dim (foutid, 'lat', nlat, latid) + if (status .ne. NF_NOERR) call handle_err(status) +!=====Jinbo Xie===== + status = nf_def_dim (foutid, 'nvar_dirOA', nvar_dirOA, varid) + if (status .ne. NF_NOERR) call handle_err(status) + status = nf_def_dim (foutid, 'nvar_dirOL', nvar_dirOL, var2id) + if (status .ne. NF_NOERR) call handle_err(status) + status = nf_def_dim (foutid, 'indexb', indexb, indexbid) + if (status .ne. NF_NOERR) call handle_err(status) +!=====Jinbo Xie===== + ! + ! Create variable for output + ! + print *,"Create variable for output" + +!========Jinbo Xie======== + ocdim(1)=lonid + ocdim(2)=latid + status = nf_def_var (foutid,'OC', NF_DOUBLE, 2, ocdim, ocid) + oadim(1)=lonid + oadim(2)=latid + oadim(3)=varid + status = nf_def_var (foutid,'OA', NF_DOUBLE, 3, oadim, oaid) + oldim(1)=lonid + oldim(2)=latid + oldim(3)=var2id + status = nf_def_var (foutid,'OL', NF_DOUBLE, 3, oldim, olid) + terroutdim(1)=lonid + terroutdim(2)=latid + terroutdim(3)=indexbid + !name + terroutchar(1)="terr" + terroutchar(2)="terrx" + terroutchar(3)="terry" + terroutchar(4)="wt" +!#if 0 + do i=1,4 + status = nf_def_var (foutid, terroutchar(i), NF_DOUBLE, 3, & + terroutdim, terroutid(i)) + enddo +!#endif + !dxy + status = nf_def_var (foutid,'dxy', NF_DOUBLE, 3, oldim, dxyid) +!#endif + +#if 0 + status = nf_def_var (foutid,'OL1', NF_DOUBLE, 2, ocdim, ol1id) + status = nf_def_var (foutid,'OL2', NF_DOUBLE, 2, ocdim, ol2id) + status = nf_def_var (foutid,'OL3', NF_DOUBLE, 2, ocdim, ol3id) + status = nf_def_var (foutid,'OL4', NF_DOUBLE, 2, ocdim, ol4id) + status = nf_def_var (foutid,'OA1', NF_DOUBLE, 2, ocdim, oa1id) + status = nf_def_var (foutid,'OA2', NF_DOUBLE, 2, ocdim, oa2id) + status = nf_def_var (foutid,'OA3', NF_DOUBLE, 2, ocdim, oa3id) + status = nf_def_var (foutid,'OA4', NF_DOUBLE, 2, ocdim, oa4id) +#endif +!========Jinbo Xie======== + + + htopodim(1)=lonid + htopodim(2)=latid + + if (lprepare_fv_smoothing_routine) then + status = nf_def_var (foutid,'htopo', NF_DOUBLE, 2, htopodim, terrid) + else + status = nf_def_var (foutid,'PHIS', NF_DOUBLE, 2, htopodim, terrid) + end if + if (status .ne. NF_NOERR) call handle_err(status) + + landfdim(1)=lonid + landfdim(2)=latid + + if (lprepare_fv_smoothing_routine) then + status = nf_def_var (foutid,'ftopo', NF_DOUBLE, 2, landfdim, landfracid) + else + status = nf_def_var (foutid,'LANDFRAC', NF_DOUBLE, 2, landfdim, landfracid) + end if + + if (status .ne. NF_NOERR) call handle_err(status) + + sghdim(1)=lonid + sghdim(2)=latid + + status = nf_def_var (foutid,'SGH', NF_DOUBLE, 2, sghdim, sghid) + if (status .ne. NF_NOERR) call handle_err(status) + + sgh30dim(1)=lonid + sgh30dim(2)=latid + + status = nf_def_var (foutid,'SGH30', NF_DOUBLE, 2, sgh30dim, sgh30id) + if (status .ne. NF_NOERR) call handle_err(status) + + landmcoslatdim(1)=lonid + landmcoslatdim(2)=latid + + status = nf_def_var (foutid,'LANDM_COSLAT', NF_DOUBLE, 2, landmcoslatdim, landm_coslatid) + if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_def_var (foutid,'lat', NF_DOUBLE, 1, latid, latvid) + if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_def_var (foutid,'lon', NF_DOUBLE, 1, lonid, lonvid) + if (status .ne. NF_NOERR) call handle_err(status) + + ! + ! Create attributes for output variables + ! + status = nf_put_att_text (foutid,terrid,'long_name', 21, 'surface geopotential') + status = nf_put_att_text (foutid,terrid,'units', 5, 'm2/s2') + status = nf_put_att_text (foutid,terrid,'filter', 35, 'area averaged from ncube3000 data') + status = nf_put_att_double (foutid, terrid, 'missing_value', nf_double, 1, fillvalue) + status = nf_put_att_double (foutid, terrid, '_FillValue' , nf_double, 1, fillvalue) + + + status = nf_put_att_double (foutid, sghid, 'missing_value', nf_double, 1, fillvalue) + status = nf_put_att_double (foutid, sghid, '_FillValue' , nf_double, 1, fillvalue) + status = nf_put_att_text (foutid, sghid, 'long_name' , 48, & + 'standard deviation of 3km cubed-sphere elevation and target grid elevation') + status = nf_put_att_text (foutid, sghid, 'units' , 1, 'm') + status = nf_put_att_text (foutid, sghid, 'filter' , 4, 'none') + + status = nf_put_att_double (foutid, sgh30id, 'missing_value', nf_double, 1, fillvalue) + status = nf_put_att_double (foutid, sgh30id, '_FillValue' , nf_double, 1, fillvalue) + status = nf_put_att_text (foutid, sgh30id, 'long_name' , 49, & + 'standard deviation of 30s elevation from 3km cubed-sphere cell average height') + status = nf_put_att_text (foutid, sgh30id, 'units' , 1, 'm') + status = nf_put_att_text (foutid, sgh30id, 'filter' , 4, 'none') + + status = nf_put_att_double (foutid, landm_coslatid, 'missing_value', nf_double, 1, fillvalue) + status = nf_put_att_double (foutid, landm_coslatid, '_FillValue' , nf_double, 1, fillvalue) + status = nf_put_att_text (foutid, landm_coslatid, 'long_name' , 23, 'smoothed land fraction') + status = nf_put_att_text (foutid, landm_coslatid, 'filter' , 4, 'none') + + status = nf_put_att_double (foutid, landfracid, 'missing_value', nf_double, 1, fillvalue) + status = nf_put_att_double (foutid, landfracid, '_FillValue' , nf_double, 1, fillvalue) + status = nf_put_att_text (foutid, landfracid, 'long_name', 21, 'gridbox land fraction') + status = nf_put_att_text (foutid, landfracid, 'filter', 40, 'area averaged from 30-sec USGS raw data') + + + status = nf_put_att_text (foutid,latvid,'long_name', 8, 'latitude') + if (status .ne. NF_NOERR) call handle_err(status) + status = nf_put_att_text (foutid,latvid,'units', 13, 'degrees_north') + if (status .ne. NF_NOERR) call handle_err(status) + ! status = nf_put_att_text (foutid,latvid,'units', 21, 'cell center locations') + ! if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_put_att_text (foutid,lonvid,'long_name', 9, 'longitude') + if (status .ne. NF_NOERR) call handle_err(status) + status = nf_put_att_text (foutid,lonvid,'units', 12, 'degrees_east') + if (status .ne. NF_NOERR) call handle_err(status) + ! status = nf_put_att_text (foutid,lonvid,'units' , 21, 'cell center locations') + ! if (status .ne. NF_NOERR) call handle_err(status) + + status = nf_put_att_text (foutid,NF_GLOBAL,'source', 27, 'USGS 30-sec dataset GTOPO30') + if (status .ne. NF_NOERR) call handle_err(status) + status = nf_put_att_text (foutid,NF_GLOBAL,'title', 24, '30-second USGS topo data') + if (status .ne. NF_NOERR) call handle_err(status) + call DATE_AND_TIME(DATE=datestring) + status = nf_put_att_text (foutid,NF_GLOBAL,'history',25, 'Written on date: ' // datestring ) + if (status .ne. NF_NOERR) call handle_err(status) + + !=====Jinbo Xie===== + status = nf_put_att_text (foutid,oaid,'note', 40, '(2)+1 in nvar_dirOA to avoid bug in io') + do i=1,4 + status = nf_put_att_double (foutid, terroutid(i),& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, terroutid(i),& + '_FillValue' , nf_double, 1,fillvalue) + enddo + + status = nf_put_att_double (foutid, oa1id,& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, oa1id,& + '_FillValue' , nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, oa2id,& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, oa2id,& + '_FillValue' , nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, oa3id,& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, oa3id,& + '_FillValue' , nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, oa4id,& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, oa4id,& + '_FillValue' , nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, ol1id,& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, ol1id,& + '_FillValue' , nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, ol2id,& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, ol2id,& + '_FillValue' , nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, ol3id,& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, ol3id,& + '_FillValue' , nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, ol4id,& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, ol4id,& + '_FillValue' , nf_double, 1,fillvalue) + !=====Jinbo Xie===== + + ! + ! End define mode for output file + ! + status = nf_enddef (foutid) + if (status .ne. NF_NOERR) call handle_err(status) + ! + ! Write variable for output +!========Jinbo Xie======== +print*,"writing oc data",MINVAL(oc),MAXVAL(oc) +status = nf_put_var_double (foutid, ocid, oc) +if (status .ne. NF_NOERR) call handle_err(status) +!oa,ol +print*,"writing oa data",MINVAL(oa),MAXVAL(oa) +status = nf_put_var_double (foutid, oaid, oa) +if (status .ne. NF_NOERR) call handle_err(status) +print*,"writing ol data",MINVAL(ol),MAXVAL(ol) +status = nf_put_var_double (foutid, olid, ol) + +!============ +#if 0 +print*,"writing oa1 data",MINVAL(oa),MAXVAL(oa) +status = nf_put_var_double (foutid, oa1id, oa(:,1)) +if (status .ne. NF_NOERR) call handle_err(status) +print*,"writing ol1 data",MINVAL(ol),MAXVAL(ol) +status = nf_put_var_double (foutid, ol1id, ol(:,1)) +print*,"writing oa2 data",MINVAL(oa),MAXVAL(oa) +status = nf_put_var_double (foutid, oa2id, oa(:,2)) +if (status .ne. NF_NOERR) call handle_err(status) +print*,"writing ol2 data",MINVAL(ol),MAXVAL(ol) +status = nf_put_var_double (foutid, ol2id, ol(:,2)) +print*,"writing oa3 data",MINVAL(oa),MAXVAL(oa) +status = nf_put_var_double (foutid, oa3id, oa(:,3)) +if (status .ne. NF_NOERR) call handle_err(status) +print*,"writing ol3 data",MINVAL(ol),MAXVAL(ol) +status = nf_put_var_double (foutid, ol3id, ol(:,3)) +print*,"writing oa4 data",MINVAL(oa),MAXVAL(oa) +status = nf_put_var_double (foutid, oa4id, oa(:,4)) +if (status .ne. NF_NOERR) call handle_err(status) +print*,"writing ol4 data",MINVAL(ol),MAXVAL(ol) +status = nf_put_var_double (foutid, ol4id, ol(:,4)) +#endif +!=========== + + +if (status .ne. NF_NOERR) call handle_err(status) +!#if 0 + do i=1,4 + status = nf_put_att_double (foutid, terroutid(i),& + 'missing_value', nf_double, 1,fillvalue) + status = nf_put_att_double (foutid, terroutid(i),& + '_FillValue' , nf_double, 1,fillvalue) + print*,"writing"//terroutchar(i)//" data",& + MINVAL(terrout(i,:,:)),MAXVAL(terrout(i,:,:)) + status = nf_put_var_double (foutid, terroutid(i), terrout(i,:,:)) + if (status .ne. NF_NOERR) call handle_err(status) + enddo +!#endif + +!#if 0 + print*,"writing dxy data",MINVAL(dxy),MAXVAL(dxy) + status = nf_put_var_double (foutid, dxyid, dxy) + if (status .ne. NF_NOERR) call handle_err(status) +!#endif +!========Jinbo Xie======== + ! + print*,"writing terrain data",MINVAL(terr),MAXVAL(terr) + if (lprepare_fv_smoothing_routine) then + status = nf_put_var_double (foutid, terrid, terr) + else + status = nf_put_var_double (foutid, terrid, terr*9.80616) + end if + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing terrain data" + + print*,"writing landfrac data",MINVAL(landfrac),MAXVAL(landfrac) + status = nf_put_var_double (foutid, landfracid, landfrac) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing landfrac data" + + print*,"writing sgh data",MINVAL(sgh),MAXVAL(sgh) + status = nf_put_var_double (foutid, sghid, sgh) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing sgh data" + + print*,"writing sgh30 data",MINVAL(sgh30),MAXVAL(sgh30) + status = nf_put_var_double (foutid, sgh30id, sgh30) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing sgh30 data" + + print*,"writing landm_coslat data",MINVAL(landm_coslat),MAXVAL(landm_coslat) + status = nf_put_var_double (foutid, landm_coslatid, landm_coslat) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing sgh30 data" + ! + print*,"writing lat data" + status = nf_put_var_double (foutid, latvid, latar) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing lat data" + + print*,"writing lon data" + status = nf_put_var_double (foutid, lonvid, lonar) + if (status .ne. NF_NOERR) call handle_err(status) + print*,"done writing lon data" + ! + ! Close output file + ! + print *,"close file" + status = nf_close (foutid) + if (status .ne. NF_NOERR) call handle_err(status) +end subroutine wrtncdf_rll +!************************************************************************ +!!handle_err +!************************************************************************ +! +!!ROUTINE: handle_err +!!DESCRIPTION: error handler +!-------------------------------------------------------------------------- + +subroutine handle_err(status) + + implicit none + +# include + + integer status + + if (status .ne. nf_noerr) then + print *, nf_strerror(status) + stop 'Stopped' + endif + +end subroutine handle_err + + +SUBROUTINE coarsen(f,fcoarse,nf,n,dA_coarse) + use shr_kind_mod, only: r8 => shr_kind_r8 + IMPLICIT NONE + REAL (R8), DIMENSION(n) , INTENT(IN) :: f + REAL (R8), DIMENSION(n/nf), INTENT(OUT) :: fcoarse + INTEGER, INTENT(in) :: n,nf + REAL(R8), DIMENSION(INT(SQRT(DBLE(n/6)))/nf,INT(SQRT(DBLE(n/6)))/nf),INTENT(OUT) :: dA_coarse + !must be an even number + ! + ! local workspace + ! + ! ncube = INT(SQRT(DBLE(n/6))) + + REAL(R8), DIMENSION(INT(SQRT(DBLE(n/6))),INT(SQRT(DBLE(n/6)))):: dA + REAL (R8) :: sum, sum_area,tmp + INTEGER :: jx,jy,jp,ii,ii_coarse,coarse_ncube,ncube + INTEGER :: jx_coarse,jy_coarse,jx_s,jy_s + + + ! REAL(R8), DIMENSION(INT(SQRT(DBLE(n/6)))/nf,INT(SQRT(DBLE(n/6)))/nf) :: dAtmp + + ncube = INT(SQRT(DBLE(n/6))) + coarse_ncube = ncube/nf + + IF (ABS(DBLE(ncube)/DBLE(nf)-coarse_ncube)>0.000001) THEN + WRITE(*,*) "ncube/nf must be an integer" + WRITE(*,*) "ncube and nf: ",ncube,nf + STOP + END IF + + da_coarse = 0.0 + + WRITE(*,*) "compute all areas" + CALL EquiangularAllAreas(ncube, dA) + ! CALL EquiangularAllAreas(coarse_ncube, dAtmp)!dbg + tmp = 0.0 + DO jp=1,6 + DO jy_coarse=1,coarse_ncube + DO jx_coarse=1,coarse_ncube + ! + ! inner loop + ! + sum = 0.0 + sum_area = 0.0 + DO jy_s=1,nf + jy = (jy_coarse-1)*nf+jy_s + DO jx_s=1,nf + jx = (jx_coarse-1)*nf+jx_s + ii = (jp-1)*ncube*ncube+(jy-1)*ncube+jx + sum = sum +f(ii)*dA(jx,jy) + sum_area = sum_area+dA(jx,jy) + ! WRITE(*,*) "jx,jy",jx,jy + END DO + END DO + tmp = tmp+sum_area + da_coarse(jx_coarse,jy_coarse) = sum_area + ! WRITE(*,*) "jx_coarse,jy_coarse",jx_coarse,jy_coarse,& + ! da_coarse(jx_coarse,jy_coarse)-datmp(jx_coarse,jy_coarse) + ii_coarse = (jp-1)*coarse_ncube*coarse_ncube+(jy_coarse-1)*coarse_ncube+jx_coarse + fcoarse(ii_coarse) = sum/sum_area + END DO + END DO + END DO + WRITE(*,*) "coarsened surface area",tmp-4.0*3.141592654 +END SUBROUTINE COARSEN + +SUBROUTINE overlap_weights(weights_lgr_index_all,weights_eul_index_all,weights_all,& + jall,ncube,ngauss,ntarget,ncorner,jmax_segments,target_corner_lon,target_corner_lat,nreconstruction) + use shr_kind_mod, only: r8 => shr_kind_r8 + use remap + IMPLICIT NONE + + + INTEGER, INTENT(INOUT) :: jall !anticipated number of weights + INTEGER, INTENT(IN) :: ncube, ngauss, ntarget, jmax_segments, ncorner, nreconstruction + + INTEGER, DIMENSION(jall,3), INTENT(OUT) :: weights_eul_index_all + REAL(R8), DIMENSION(jall,nreconstruction) , INTENT(OUT) :: weights_all + INTEGER, DIMENSION(jall) , INTENT(OUT) :: weights_lgr_index_all + + REAL(R8), DIMENSION(ncorner,ntarget), INTENT(IN) :: target_corner_lon, target_corner_lat + + INTEGER, DIMENSION(ncorner+1) :: ipanel_array, ipanel_tmp + REAL(R8), DIMENSION(ncorner) :: lat, lon + REAL(R8), DIMENSION(0:ncube+2):: xgno, ygno + REAL(R8), DIMENSION(0:ncorner+1) :: xcell, ycell + + REAL(R8), DIMENSION(ngauss) :: gauss_weights, abscissae + + REAL(R8) :: da, tmp, alpha, beta + REAL (r8), PARAMETER :: pi = 3.14159265358979323846264338327 + REAL (r8), PARAMETER :: piq = 0.25*pi + REAL (r8), PARAMETER :: pih = 0.50*pi + INTEGER :: i, j,ncorner_this_cell,k,ip,ipanel,ii,jx,jy,jcollect + integer :: alloc_error + + REAL (r8), PARAMETER :: rad2deg = 180.0/pi + + real(r8), allocatable, dimension(:,:) :: weights + integer , allocatable, dimension(:,:) :: weights_eul_index + + + LOGICAL:: ldbg = .FAlSE. + + INTEGER :: jall_anticipated + + jall_anticipated = jall + + ipanel_array = -99 + ! + da = pih/DBLE(ncube) + xgno(0) = -bignum + DO i=1,ncube+1 + xgno(i) = TAN(-piq+(i-1)*da) + END DO + xgno(ncube+2) = bignum + ygno = xgno + + CALL glwp(ngauss,gauss_weights,abscissae) + + + allocate (weights(jmax_segments,nreconstruction),stat=alloc_error ) + allocate (weights_eul_index(jmax_segments,2),stat=alloc_error ) + + tmp = 0.0 + jall = 1 + DO i=1,ntarget + WRITE(*,*) "cell",i," ",100.0*DBLE(i)/DBLE(ntarget),"% done" + ! + !--------------------------------------------------- + ! + ! determine how many vertices the cell has + ! + !--------------------------------------------------- + ! + CALL remove_duplicates_latlon(ncorner,target_corner_lon(:,i),target_corner_lat(:,i),& + ncorner_this_cell,lon,lat,1.0E-10,ldbg) + + IF (ldbg) THEN + WRITE(*,*) "number of vertices ",ncorner_this_cell + WRITE(*,*) "vertices locations lon,",lon(1:ncorner_this_cell)*rad2deg + WRITE(*,*) "vertices locations lat,",lat(1:ncorner_this_cell)*rad2deg + DO j=1,ncorner_this_cell + WRITE(*,*) lon(j)*rad2deg, lat(j)*rad2deg + END DO + WRITE(*,*) " " + END IF + ! + !--------------------------------------------------- + ! + ! determine how many and which panels the cell spans + ! + !--------------------------------------------------- + ! + DO j=1,ncorner_this_cell + CALL CubedSphereABPFromRLL(lon(j), lat(j), alpha, beta, ipanel_tmp(j), .TRUE.) + IF (ldbg) WRITE(*,*) "ipanel for corner ",j," is ",ipanel_tmp(j) + END DO + ipanel_tmp(ncorner_this_cell+1) = ipanel_tmp(1) + ! make sure to include possible overlap areas not on the face the vertices are located + IF (MINVAL(lat(1:ncorner_this_cell))<-pi/6.0) THEN + ! include South-pole panel in search + ipanel_tmp(ncorner_this_cell+1) = 5 + IF (ldbg) WRITE(*,*) "add panel 5 to search" + END IF + IF (MAXVAL(lat(1:ncorner_this_cell))>pi/6.0) THEN + ! include North-pole panel in search + ipanel_tmp(ncorner_this_cell+1) = 6 + IF (ldbg) WRITE(*,*) "add panel 6 to search" + END IF + ! + ! remove duplicates in ipanel_tmp + ! + CALL remove_duplicates_integer(ncorner_this_cell+1,ipanel_tmp(1:ncorner_this_cell+1),& + k,ipanel_array(1:ncorner_this_cell+1)) + ! + !--------------------------------------------------- + ! + ! loop over panels with possible overlap areas + ! + !--------------------------------------------------- + ! + DO ip = 1,k + ipanel = ipanel_array(ip) + DO j=1,ncorner_this_cell + ii = ipanel + CALL CubedSphereABPFromRLL(lon(j), lat(j), alpha, beta, ii,.FALSE.) + IF (j==1) THEN + jx = CEILING((alpha + piq) / da) + jy = CEILING((beta + piq) / da) + END IF + xcell(ncorner_this_cell+1-j) = TAN(alpha) + ycell(ncorner_this_cell+1-j) = TAN(beta) + END DO + xcell(0) = xcell(ncorner_this_cell) + ycell(0) = ycell(ncorner_this_cell) + xcell(ncorner_this_cell+1) = xcell(1) + ycell(ncorner_this_cell+1) = ycell(1) + + jx = MAX(MIN(jx,ncube+1),0) + jy = MAX(MIN(jy,ncube+1),0) + + CALL compute_weights_cell(xcell(0:ncorner_this_cell+1),ycell(0:ncorner_this_cell+1),& + jx,jy,nreconstruction,xgno,ygno,& + 1, ncube+1, 1,ncube+1, tmp,& + ngauss,gauss_weights,abscissae,weights,weights_eul_index,jcollect,jmax_segments,& + ncube,0,ncorner_this_cell,ldbg) + + weights_all(jall:jall+jcollect-1,1:nreconstruction) = weights(1:jcollect,1:nreconstruction) + + weights_eul_index_all(jall:jall+jcollect-1,1:2) = weights_eul_index(1:jcollect,:) + weights_eul_index_all(jall:jall+jcollect-1, 3) = ipanel + weights_lgr_index_all(jall:jall+jcollect-1 ) = i + + jall = jall+jcollect + IF (jall>jall_anticipated) THEN + WRITE(*,*) "more weights than anticipated" + WRITE(*,*) "increase jall" + STOP + END IF + IF (ldbg) WRITE(*,*) "jcollect",jcollect + END DO + END DO + jall = jall-1 + WRITE(*,*) "sum of all weights divided by surface area of sphere =",tmp/(4.0*pi) + WRITE(*,*) "actual number of weights",jall + WRITE(*,*) "anticipated number of weights",jall_anticipated + IF (jall>jall_anticipated) THEN + WRITE(*,*) "anticipated number of weights < actual number of weights" + WRITE(*,*) "increase jall!" + STOP + END IF + WRITE(*,*) MINVAL(weights_all(1:jall,1)),MAXVAL(weights_all(1:jall,1)) + IF (ABS(tmp/(4.0*pi))-1.0>0.001) THEN + WRITE(*,*) "sum of all weights does not match the surface area of the sphere" + WRITE(*,*) "sum of all weights is : ",tmp + WRITE(*,*) "surface area of sphere: ",4.0*pi + STOP + END IF +END SUBROUTINE overlap_weights + + +!------------------------------------------------------------------------------ +! SUBROUTINE CubedSphereABPFromRLL +! +! Description: +! Determine the (alpha,beta,panel) coordinate of a point on the sphere from +! a given regular lat lon coordinate. +! +! Parameters: +! lon - Coordinate longitude +! lat - Coordinate latitude +! alpha (OUT) - Alpha coordinate +! beta (OUT) - Beta coordinate +! ipanel (OUT) - Face panel +!------------------------------------------------------------------------------ +SUBROUTINE CubedSphereABPFromRLL(lon, lat, alpha, beta, ipanel, ldetermine_panel) + use shr_kind_mod, only: r8 => shr_kind_r8 + IMPLICIT NONE + + REAL (R8), INTENT(IN) :: lon, lat + REAL (R8), INTENT(OUT) :: alpha, beta + INTEGER :: ipanel + LOGICAL, INTENT(IN) :: ldetermine_panel + REAL (r8), PARAMETER :: pi = 3.14159265358979323846264338327 + REAL (r8), PARAMETER :: piq = 0.25*pi + REAL (r8), PARAMETER :: rotate_cube = 0.0 + + ! Local variables + REAL (R8) :: xx, yy, zz, pm + REAL (R8) :: sx, sy, sz + INTEGER :: ix, iy, iz + + ! Translate to (x,y,z) space + xx = COS(lon-rotate_cube) * COS(lat) + yy = SIN(lon-rotate_cube) * COS(lat) + zz = SIN(lat) + + pm = MAX(ABS(xx), ABS(yy), ABS(zz)) + + ! Check maximality of the x coordinate + IF (pm == ABS(xx)) THEN + IF (xx > 0) THEN; ix = 1; ELSE; ix = -1; ENDIF + ELSE + ix = 0 + ENDIF + + ! Check maximality of the y coordinate + IF (pm == ABS(yy)) THEN + IF (yy > 0) THEN; iy = 1; ELSE; iy = -1; ENDIF + ELSE + iy = 0 + ENDIF + + ! Check maximality of the z coordinate + IF (pm == ABS(zz)) THEN + IF (zz > 0) THEN; iz = 1; ELSE; iz = -1; ENDIF + ELSE + iz = 0 + ENDIF + + ! Panel assignments + IF (ldetermine_panel) THEN + IF (iz == 1) THEN + ipanel = 6; sx = yy; sy = -xx; sz = zz + + ELSEIF (iz == -1) THEN + ipanel = 5; sx = yy; sy = xx; sz = -zz + + ELSEIF ((ix == 1) .AND. (iy /= 1)) THEN + ipanel = 1; sx = yy; sy = zz; sz = xx + + ELSEIF ((ix == -1) .AND. (iy /= -1)) THEN + ipanel = 3; sx = -yy; sy = zz; sz = -xx + + ELSEIF ((iy == 1) .AND. (ix /= -1)) THEN + ipanel = 2; sx = -xx; sy = zz; sz = yy + + ELSEIF ((iy == -1) .AND. (ix /= 1)) THEN + ipanel = 4; sx = xx; sy = zz; sz = -yy + + ELSE + WRITE(*,*) 'Fatal Error: CubedSphereABPFromRLL failed' + WRITE(*,*) '(xx, yy, zz) = (', xx, ',', yy, ',', zz, ')' + WRITE(*,*) 'pm =', pm, ' (ix, iy, iz) = (', ix, ',', iy, ',', iz, ')' + STOP + ENDIF + ELSE + IF (ipanel == 6) THEN + sx = yy; sy = -xx; sz = zz + ELSEIF (ipanel == 5) THEN + sx = yy; sy = xx; sz = -zz + ELSEIF (ipanel == 1) THEN + sx = yy; sy = zz; sz = xx + ELSEIF (ipanel == 3) THEN + sx = -yy; sy = zz; sz = -xx + ELSEIF (ipanel == 2) THEN + sx = -xx; sy = zz; sz = yy + ELSEIF (ipanel == 4) THEN + sx = xx; sy = zz; sz = -yy + ELSE + WRITE(*,*) "ipanel out of range",ipanel + STOP + END IF + END IF + + ! Use panel information to calculate (alpha, beta) coords + alpha = ATAN(sx / sz) + beta = ATAN(sy / sz) + +END SUBROUTINE CubedSphereABPFromRLL + +!------------------------------------------------------------------------------ +! SUBROUTINE EquiangularAllAreas +! +! Description: +! Compute the area of all cubed sphere grid cells, storing the results in +! a two dimensional array. +! +! Parameters: +! icube - Resolution of the cubed sphere +! dA (OUT) - Output array containing the area of all cubed sphere grid cells +!------------------------------------------------------------------------------ +SUBROUTINE EquiangularAllAreas(icube, dA) + use shr_kind_mod, only: r8 => shr_kind_r8 + IMPLICIT NONE + + INTEGER, INTENT(IN) :: icube + REAL (r8), DIMENSION(icube,icube), INTENT(OUT) :: dA + + ! Local variables + INTEGER :: k, k1, k2 + REAL (r8) :: a1, a2, a3, a4 + REAL (r8), DIMENSION(icube+1,icube+1) :: ang + REAL (r8), DIMENSION(icube+1) :: gp + + REAL (r8), PARAMETER :: pi = 3.14159265358979323846264338327 + REAL (r8), PARAMETER :: piq = 0.25*pi + + + !#ifdef DBG + REAL (r8) :: dbg1 !DBG + !#endif + + ! Recall that we are using equi-angular spherical gridding + ! Compute the angle between equiangular cubed sphere projection grid lines. + DO k = 1, icube+1 + gp(k) = -piq + (pi/DBLE(2*(icube))) * DBLE(k-1) + ENDDO + + DO k2=1,icube+1 + DO k1=1,icube+1 + ang(k1,k2) =ACOS(-SIN(gp(k1)) * SIN(gp(k2))) + ENDDO + ENDDO + + DO k2=1,icube + DO k1=1,icube + a1 = ang(k1 , k2 ) + a2 = pi - ang(k1+1, k2 ) + a3 = pi - ang(k1 , k2+1) + a4 = ang(k1+1, k2+1) + ! area = r*r*(-2*pi+sum(interior angles)) + DA(k1,k2) = -2.0*pi+a1+a2+a3+a4 + ENDDO + ENDDO + + !#ifdef DBG + ! Only for debugging - test consistency + dbg1 = 0.0 !DBG + DO k2=1,icube + DO k1=1,icube + dbg1 = dbg1 + DA(k1,k2) !DBG + ENDDO + ENDDO + write(*,*) 'DAcube consistency: ',dbg1-4.0*pi/6.0 !DBG + !#endif +END SUBROUTINE EquiangularAllAreas + + +!------------------------------------------------------------------------------ +! SUBROUTINE CubedSphereRLLFromABP +! +! Description: +! Determine the lat lon coordinate of a point on a sphere given its +! (alpha,beta,panel) coordinate. +! +! Parameters: +! alpha - Alpha coordinate +! beta - Beta coordinate +! panel - Cubed sphere panel id +! lon (OUT) - Calculated longitude +! lat (OUT) - Calculated latitude +!------------------------------------------------------------------------------ +SUBROUTINE CubedSphereRLLFromABP(alpha, beta, ipanel, lon, lat) + use shr_kind_mod, only: r8 => shr_kind_r8 + IMPLICIT NONE + REAL (r8), INTENT(IN) :: alpha, beta + INTEGER , INTENT(IN) :: ipanel + REAL (r8), INTENT(OUT) :: lon, lat + ! Local variables + REAL (r8) :: xx, yy, zz, rotate_cube + REAL (r8), PARAMETER :: pi = 3.14159265358979323846264338327 + REAL (r8), PARAMETER :: piq = 0.25*pi + + rotate_cube = 0.0 + ! Convert to cartesian coordinates + CALL CubedSphereXYZFromABP(alpha, beta, ipanel, xx, yy, zz) + ! Convert back to lat lon + lat = ASIN(zz) + if (xx==0.0.and.yy==0.0) THEN + lon = 0.0 + else + lon = ATAN2(yy, xx) +rotate_cube + IF (lon<0.0) lon=lon+2.0*pi + IF (lon>2.0*pi) lon=lon-2.0*pi + end if +END SUBROUTINE CubedSphereRLLFromABP + +!------------------------------------------------------------------------------ +! SUBROUTINE CubedSphereXYZFromABP +! +! Description: +! Determine the Cartesian coordinate of a point on a sphere given its +! (alpha,beta,panel) coordinate. +! +! Parameters: +! alpha - Alpha coordinate +! beta - Beta coordinate +! panel - Cubed sphere panel id +! xx (OUT) - Calculated x coordinate +! yy (OUT) - Calculated y coordinate +! zz (OUT) - Calculated z coordinate +!------------------------------------------------------------------------------ +SUBROUTINE CubedSphereXYZFromABP(alpha, beta, ipanel, xx, yy, zz) + use shr_kind_mod, only: r8 => shr_kind_r8 + IMPLICIT NONE + + REAL (r8), INTENT(IN) :: alpha, beta + INTEGER , INTENT(IN) :: ipanel + REAL (r8), INTENT(OUT) :: xx, yy, zz + ! Local variables + REAL (r8) :: a1, b1, pm + REAL (r8) :: sx, sy, sz + + ! Convert to Cartesian coordinates + a1 = TAN(alpha) + b1 = TAN(beta) + + sz = (1.0 + a1 * a1 + b1 * b1)**(-0.5) + sx = sz * a1 + sy = sz * b1 + ! Panel assignments + IF (ipanel == 6) THEN + yy = sx; xx = -sy; zz = sz + ELSEIF (ipanel == 5) THEN + yy = sx; xx = sy; zz = -sz + ELSEIF (ipanel == 1) THEN + yy = sx; zz = sy; xx = sz + ELSEIF (ipanel == 3) THEN + yy = -sx; zz = sy; xx = -sz + ELSEIF (ipanel == 2) THEN + xx = -sx; zz = sy; yy = sz + ELSEIF (ipanel == 4) THEN + xx = sx; zz = sy; yy = -sz + ELSE + WRITE(*,*) 'Fatal Error: Panel out of range in CubedSphereXYZFromABP' + WRITE(*,*) '(alpha, beta, panel) = (', alpha, ',', beta, ',', ipanel, ')' + STOP + ENDIF +END SUBROUTINE CubedSphereXYZFromABP + + +SUBROUTINE remove_duplicates_integer(n_in,f_in,n_out,f_out) + use shr_kind_mod, only: r8 => shr_kind_r8 + integer, intent(in) :: n_in + integer,dimension(n_in), intent(in) :: f_in + integer, intent(out) :: n_out + integer,dimension(n_in), intent(out) :: f_out + ! + ! local work space + ! + integer :: k,i,j + ! + ! remove duplicates in ipanel_tmp + ! + k = 1 + f_out(1) = f_in(1) + outer: do i=2,n_in + do j=1,k + ! if (f_out(j) == f_in(i)) then + if (ABS(f_out(j)-f_in(i))<1.0E-10) then + ! Found a match so start looking again + cycle outer + end if + end do + ! No match found so add it to the output + k = k + 1 + f_out(k) = f_in(i) + end do outer + n_out = k +END SUBROUTINE remove_duplicates_integer + +SUBROUTINE remove_duplicates_latlon(n_in,lon_in,lat_in,n_out,lon_out,lat_out,tiny,ldbg) + use shr_kind_mod, only: r8 => shr_kind_r8 + integer, intent(in) :: n_in + real(r8),dimension(n_in), intent(inout) :: lon_in,lat_in + real, intent(in) :: tiny + integer, intent(out) :: n_out + real(r8),dimension(n_in), intent(out) :: lon_out,lat_out + logical :: ldbg + ! + ! local work space + ! + integer :: k,i,j + REAL (r8), PARAMETER :: pi = 3.14159265358979323846264338327 + REAL (r8), PARAMETER :: pih = 0.50*pi + ! + ! for pole points: make sure the longitudes are identical so that algorithm below works properly + ! + do i=2,n_in + if (abs(lat_in(i)-pih) ' + print *, ' ' + print *, 'REQUIRED ARGUMENTS: ' + print *, ' --target-grid Target grid descriptor in SCRIP format ' + print *, ' --input-topography Input USGS topography on cube sphere ' + print *, ' --output-topography Output topography on target grid ' + print *, ' ' + print *, 'OPTIONAL ARGUMENTS: ' + print *, ' --smoothed-topography Input smoothed topography (for surface ' + print *, ' roughness calculation). If present, ' + print *, ' output will contain estimate of subgrid' + print *, ' surface roughness (SGH). Note that the ' + print *, ' variance in elevation from the 30s to ' + print *, ' 3km grid (SGH30) is also downscaled, ' + print *, ' but does not depend on the smoothing. ' + print *, ' ' + print *, 'DESCRIPTION: ' + print *, 'This code performs rigorous remapping of topography variables on a cubed- ' + print *, 'sphere grid to any target grid. The code is documented in: ' + print *, ' ' + print *, ' Lauritzen, Nair and Ullrich, 2010, J. Comput. Phys. ' + print *, ' ' + print *, 'AUTHOR: ' + print *, ' Peter Hjort Lauritzen (pel@ucar.edu), AMP/CGD/NESL/NCAR ' + print *, ' ' +end subroutine usage diff --git a/components/eam/tools/orographic_drag_toolkit/make.ncl b/components/eam/tools/orographic_drag_toolkit/make.ncl new file mode 100755 index 000000000000..d79fc234bebf --- /dev/null +++ b/components/eam/tools/orographic_drag_toolkit/make.ncl @@ -0,0 +1,21 @@ +load "/lcrc/group/e3sm/ac.xie7/Analysis/NCLep/self.ncl" +begin +vars=(/"PHIS","SGH","SGH30","LANDFRAC","LANDM_COSLAT"/) +;; +fil1="final-180-ne30pg2-mod-v3.nc" +;fil2="USGS-gtopo30_ne30np4pg2_16xdel2.c20200108.nc" +;fil3="final-180-ne30pg2.nc" +fil2="USGS-gtopo30_ne30np4pg2_x6t-SGH.c20210614.nc" +fil3="final-180-ne30pg2-v3.nc" +system("rm -r "+fil1) +system("cp -r "+fil3+" "+fil1) +;; +ff1=addfile(fil1,"w") +ff2=addfile(fil2,"r") +;; +do i=0,4 +ff1->$vars(i)$=ff2->$vars(i)$ +end do + + +end diff --git a/components/eam/tools/orographic_drag_toolkit/ogwd_sub.F90 b/components/eam/tools/orographic_drag_toolkit/ogwd_sub.F90 new file mode 100755 index 000000000000..dc3d3ba61df6 --- /dev/null +++ b/components/eam/tools/orographic_drag_toolkit/ogwd_sub.F90 @@ -0,0 +1,908 @@ +Module ogwd_sub +use shr_kind_mod, only: r8 => shr_kind_r8 +!use transform + +contains +!#if 0 +subroutine OAdir(terr,ntarget,ncube,n,nvar_dir,jall,weights_lgr_index_all,weights_eul_index_all1,weights_eul_index_all2,weights_eul_index_all3,weights_all,landfrac_target,lon_cen,lat_cen,lon_terr,lat_terr,area_target,oa_target) +!use shr_kind_mod, only: r8 => shr_kind_r8 +IMPLICIT NONE +integer ,intent(in) :: ncube,ntarget,n,nvar_dir,jall,weights_lgr_index_all(jall) +integer ,intent(in) :: weights_eul_index_all1(jall),& + weights_eul_index_all2(jall),& + weights_eul_index_all3(jall) +real(r8),intent(in) :: weights_all(jall,1),landfrac_target(ntarget) +real(r8),intent(in) :: terr(n) +!real(r8),intent(in) :: lon_cen(ntarget),& +real(r8),intent(inout) :: lon_cen(ntarget),& + lat_cen(ntarget),& + area_target(ntarget) +real(r8),intent(in) :: lon_terr(n),lat_terr(n) +real(r8),intent(out) :: oa_target(ntarget,nvar_dir) +!local +integer :: count,i,ix,iy,ip,ii,ip2,ip3 +real(r8) :: xxterr,yyterr,zzterr,ix2,iy2,xx3,yy3,zz3,ix3,iy3 +real(r8) :: wt,xhds(ntarget),yhds(ntarget),zhds(ntarget),hds(ntarget),OAx_var(ntarget),OAy_var(ntarget),OAz_var(ntarget),OA_var(ntarget) +real(r8) :: xbar(ntarget),ybar(ntarget),zbar(ntarget),lon_bar(ntarget),lat_bar(ntarget) +real(r8) :: rad,theta1 +real(r8) :: OAlon(ntarget),OAlat(ntarget),OArad(ntarget),OAx1,OAy1,OAz1 + +real(r8) :: terr_anom(n),terr_avg(ntarget),weights_ano(jall),area_target_ano(ntarget) + +xhds=0.0_r8 +yhds=0.0_r8 +zhds=0.0_r8 +hds=0.0_r8 + +xbar=0.0_r8 +ybar=0.0_r8 +zbar=0.0_r8 +lon_bar=0.0_r8 +lat_bar=0.0_r8 +OA_var=0.0_r8 +OAx_var=0.0_r8 +OAy_var=0.0_r8 +OAz_var=0.0_r8 + + +!#if 0 +terr_anom=0.0_r8 +terr_avg=0.0_r8 +do count=1,jall + i = weights_lgr_index_all(count) + ix = weights_eul_index_all1(count)!,1) + iy = weights_eul_index_all2(count)!,2) + ip = weights_eul_index_all3(count) + ! convert to 1D indexing of cubed-sphere + ii = (ip-1)*ncube*ncube+(iy-1)*ncube+ix! + wt = weights_all(count,1) + ! + terr_avg(i)=terr_avg(i)+(wt/area_target(i))*terr(ii) + !terr(ii)*wt!(wt/area_target(i))*terr(ii) +enddo + +do count=1,jall + i = weights_lgr_index_all(count) + ix = weights_eul_index_all1(count)!,1) + iy = weights_eul_index_all2(count)!,2) + ip = weights_eul_index_all3(count) + ii = (ip-1)*ncube*ncube+(iy-1)*ncube+ix + terr_anom(ii)=terr(ii)-terr_avg(i) +! +enddo +where(terr_anom.le.0) terr_anom=0.0_r8 + +do count=1,jall + i = weights_lgr_index_all(count) + ix = weights_eul_index_all1(count)!,1) + iy = weights_eul_index_all2(count)!,2) + ip = weights_eul_index_all3(count)!,3) + ! + ! convert to 1D indexing of cubed-sphere + ii = (ip-1)*ncube*ncube+(iy-1)*ncube+ix! + wt = weights_all(count,1) + rad=4.0_r8*atan(1.0_r8)/180.0_r8 + call CubedSphereABPFromRLL(lon_terr(ii)*rad,lat_terr(ii)*rad,ix2,iy2,ip2,.true.) + call CubedSphereXYZFromABP(ix2,iy2,ip2,xxterr,yyterr,zzterr) +!#if 0 + xhds(i)=xhds(i)+xxterr*terr_anom(ii)*wt + yhds(i)=yhds(i)+yyterr*terr_anom(ii)*wt + zhds(i)=zhds(i)+zzterr*terr_anom(ii)*wt + hds(i) =hds(i)+terr_anom(ii)*wt + + !masscenter for every coarse grid + !on Cartesian coord + !looking the h as ro + xbar(i)=xhds(i)/hds(i) + ybar(i)=yhds(i)/hds(i) + zbar(i)=zhds(i)/hds(i) + + call CubedSphereABPFromRLL(lon_cen(i)*rad,lat_cen(i)*rad,& + ix3,iy3,ip3,.true.) + call CubedSphereXYZFromABP(ix3,iy3,ip3,xx3,yy3,zz3) + !under Cartesian, the variability of the scale in the wind + !direction is the sqrt(x^2+y^2+z^2), the scale of the orthogonal + !3 directions + !then it is only a matter of using the original formula + !in the single direction + OA_var(i)=OA_var(i)+wt/area_target(i)& + *((xxterr-xx3)**2+(yyterr-yy3)**2+(zzterr-zz3)**2) + OAx_var(i)=OAx_var(i)+(wt/area_target(i))*(xxterr-xx3)**2 + OAy_var(i)=OAy_var(i)+(wt/area_target(i))*(yyterr-yy3)**2 + OAz_var(i)=OAz_var(i)+(wt/area_target(i))*(zzterr-zz3)**2 + OAx1=(xx3-xbar(i))/sqrt(OAx_var(i))!OA_var(i)) + OAy1=(yy3-ybar(i))/sqrt(OAy_var(i))!OA_var(i)) + OAz1=(zz3-zbar(i))/sqrt(OAz_var(i))!OA_var(i)) + !assuming a small change in lon_cen to lon_bar + !so it does not matter whether lon_cen or lon_bar + !thus we change onto lat-lon grid vector in target gridcell +#if 0 + OArad(i)= OAx1*sin(lat_cen(i)*rad)*cos(lon_cen(i)*rad)& + +OAy1*sin(lat_cen(i)*rad)*sin(lon_cen(i)*rad)& + +OAz1*cos(lat_cen(i)*rad) + OAlat(i)= OAx1*cos(lat_cen(i)*rad)*cos(lon_cen(i)*rad)& + +OAy1*cos(lat_cen(i)*rad)*sin(lon_cen(i)*rad)& + -OAz1*sin(lat_cen(i)*rad) + OAlon(i)=-OAx1*sin(lon_cen(i)*rad)& + +OAy1*cos(lon_cen(i)*rad) +#endif + !all lat_cen must use (90-lat_cen) since we only have + !latitude rather than colatitude + !this is equivalent to using induction formula sin(90-lat)=cos(lat) + !latitude is opposite of colatitude, thus OAlat is negative + OAlat(i)=-(OAx1*sin(lat_cen(i)*rad)*cos(lon_cen(i)*rad)& + +OAy1*sin(lat_cen(i)*rad)*sin(lon_cen(i)*rad)& + -OAz1*cos(lat_cen(i)*rad)) + OAlon(i)= -OAx1*sin(lon_cen(i)*rad)& + +OAy1*cos(lon_cen(i)*rad) +#if 0 + theta1=0. + oa_target(i,1) = OAlon(i)*cos(theta1*rad)+OAlat(i)*sin(theta1*rad) + theta1=90. + oa_target(i,2) = OAlon(i)*cos(theta1*rad)+OAlat(i)*sin(theta1*rad) + theta1=45. + oa_target(i,3)= OAlon(i)*cos(theta1*rad)+OAlat(i)*sin(theta1*rad) + theta1=360.-45. + oa_target(i,4)= OAlon(i)*cos(theta1*rad)+OAlat(i)*sin(theta1*rad) +#endif +!#if 0 + !reverse in order to be + !(2,ntarget),OAx,OAy + oa_target(i,1) = OAlon(i) + oa_target(i,2) = OAlat(i) + +!#endif + !landfrac may cause coastal area par to diminish + !oa_target(i,:) = oa_target(i,:)*landfrac_target(i) +enddo + !takeout abnormal values +!#if 0 + where(abs(oa_target)<.001_r8.or.& + abs(oa_target).gt.1e+7) oa_target=0.0_r8 + !where(abs(oa_target).gt.1) oa_target=1.0_r8 + where(oa_target.ne.oa_target) oa_target=0.0_r8 + +!#endif +end subroutine OAdir +!============================================================ +subroutine OAorig(terr,ntarget,ncube,n,jall,weights_lgr_index_all,weights_eul_index_all1,weights_eul_index_all2,weights_eul_index_all3,weights_all,landfrac_target,lon_terr,lat_terr,area_target,oa_target) +!use shr_kind_mod, only: r8 => shr_kind_r8 +IMPLICIT NONE +integer ,intent(in) :: ncube,ntarget,n,jall,weights_lgr_index_all(jall),weights_eul_index_all1(jall),weights_eul_index_all2(jall),weights_eul_index_all3(jall) +real(r8),intent(in) :: weights_all(jall,1),terr(n) +real(r8),intent(in) :: landfrac_target(ntarget),lon_terr(n),lat_terr(n),area_target(ntarget) +real(r8),intent(out) :: oa_target(ntarget,4) +!local +real(r8) :: xh(ntarget),yh(ntarget),height(ntarget),modexcoords(ntarget),modeycoords(ntarget),avgx(ntarget),avgy(ntarget),varx(ntarget),vary(ntarget),OAx,OAy,theta1,rad +integer(r8) :: count,i,ix,iy,ip,ii +real(r8) :: wt + + xh=0.0_r8 + yh=0.0_r8 + height=0.0_r8 + modexcoords=0.0_r8 + modeycoords=0.0_r8 + avgx=0.0_r8 + avgy=0.0_r8 + varx=0.0_r8 + vary=0.0_r8 + OAx=0.0_r8 + OAy=0.0_r8 + theta1=0.0_r8 + rad=0.0_r8 + +do count=1,jall + i = weights_lgr_index_all(count) + ix = weights_eul_index_all1(count)!,1) + iy = weights_eul_index_all2(count)!,2) + ip = weights_eul_index_all3(count)!,3) + ! + ! convert to 1D indexing of cubed-sphere + ! + ii = (ip-1)*ncube*ncube+(iy-1)*ncube+ix! + wt = weights_all(count,1) + !for OA + avgx(i)=avgx(i)+wt/area_target(i)*lon_terr(ii) + avgy(i)=avgy(i)+wt/area_target(i)*lat_terr(ii) +enddo + + +do count=1,jall + i = weights_lgr_index_all(count) + ix = weights_eul_index_all1(count)!,1) + iy = weights_eul_index_all2(count)!,2) + ip = weights_eul_index_all3(count)!,3) + ! + ! convert to 1D indexing of cubed-sphere + ! + ii = (ip-1)*ncube*ncube+(iy-1)*ncube+ix! + wt = weights_all(count,1) + !mode for both dim + xh(i)=xh(i)+wt/area_target(i)*lon_terr(ii)*terr(ii) + yh(i)=yh(i)+wt/area_target(i)*lat_terr(ii)*terr(ii) + height(i)=height(i)+wt/area_target(i)*terr(ii) + modexcoords(i)=xh(i)/(height(i))!+1e-14) + modeycoords(i)=yh(i)/(height(i))!+1e-14) + + varx(i)=varx(i)+(wt/area_target(i))*(lon_terr(ii)-avgx(i))**2 + vary(i)=vary(i)+(wt/area_target(i))*(lat_terr(ii)-avgy(i))**2 + !OAx,OAy + OAx=landfrac_target(i)*(avgx(i)-modexcoords(i))/sqrt(varx(i)) + OAy=landfrac_target(i)*(avgy(i)-modeycoords(i))/sqrt(vary(i)) + + rad=4.0*atan(1.0)/180.0 + theta1=0. + oa_target(i,1) = OAx*cos(theta1*rad)+OAy*sin(theta1*rad) + theta1=90. + oa_target(i,2) = OAx*cos(theta1*rad)+OAy*sin(theta1*rad) + theta1=45. + oa_target(i,3)= OAx*cos(theta1*rad)+OAy*sin(theta1*rad) + theta1=360.-45. + oa_target(i,4)= OAx*cos(theta1*rad)+OAy*sin(theta1*rad) + oa_target(i,:)= oa_target(i,:)*landfrac_target(i) +enddo + !takeout abnormal values + where(abs(oa_target)<.001_r8.or.abs(oa_target).gt.1e+7) oa_target=0.0 + where(abs(oa_target).gt.1) oa_target=0.0 + where(oa_target.ne.oa_target) oa_target=0.0 +end subroutine OAorig +!#endif +!=================================== +subroutine OC(terr,ntarget,ncube,n,jall,weights_lgr_index_all,weights_eul_index_all1,weights_eul_index_all2,weights_eul_index_all3,weights_all,landfrac_target,area_target,sgh_target,terr_target,oc_target) +!use shr_kind_mod, only: r8 => shr_kind_r8 +IMPLICIT NONE +integer ,intent(in) :: ncube,ntarget,n,jall,weights_lgr_index_all(jall),weights_eul_index_all1(jall),weights_eul_index_all2(jall),weights_eul_index_all3(jall) +real(r8),intent(in) :: weights_all(jall,1) +real(r8),intent(in) :: landfrac_target(ntarget),area_target(ntarget),sgh_target(ntarget),terr_target(ntarget),terr(n) +real(r8),intent(out) :: oc_target(ntarget) +!local +integer :: count,i,ix,iy,ip,ii +real(r8) :: wt + + oc_target=0.0_r8 + do count=1,jall + i = weights_lgr_index_all(count) + ix = weights_eul_index_all1(count)!,1) + iy = weights_eul_index_all2(count)!,2) + ip = weights_eul_index_all3(count)!,3) + ! convert to 1D indexing of cubed-sphere + ii = (ip-1)*ncube*ncube+(iy-1)*ncube+ix! + wt = weights_all(count,1) + oc_target(i) = oc_target(i)+(wt/area_target(i))*((terr_target(i)-terr(ii))**4)/(sgh_target(i)**4) + oc_target(i) = oc_target(i) * landfrac_target(i) + enddo + + where(abs(oc_target)<.001_r8.or.abs(oc_target).gt.1e+7) oc_target=0.0_r8 + where(abs(sgh_target).eq.0.0_r8) oc_target=0.0_r8 + where(oc_target<0.0_r8) oc_target=0.0_r8 +end subroutine OC +!======================== +subroutine OLorig(terr,ntarget,ncube,n,jall,weights_lgr_index_all,weights_eul_index_all1,weights_eul_index_all2,weights_eul_index_all3,weights_all,landfrac_target,lon_terr,lat_terr,area_target,sgh_target,target_center_lat,target_center_lon,target_corner_lat_deg,target_corner_lon_deg,ol_target) +!use shr_kind_mod, only: r8 => shr_kind_r8 +IMPLICIT NONE +integer,intent(in) :: ncube,ntarget,n,jall,weights_lgr_index_all(jall),weights_eul_index_all1(jall),weights_eul_index_all2(jall),weights_eul_index_all3(jall) +real(r8),intent(in) :: weights_all(jall,1) +real(r8),intent(in) :: landfrac_target(ntarget),area_target(ntarget),sgh_target(ntarget),terr(n),lon_terr(n),lat_terr(n) +real(r8),intent(in) :: target_center_lat(ntarget),target_center_lon(ntarget),target_corner_lat_deg(4,ntarget),target_corner_lon_deg(4,ntarget) +real(r8),intent(out) :: ol_target(ntarget,4) +!local +integer :: count,i,ix,iy,ip,ii +real(r8) :: wt,terr_if,Nw(4,ntarget),area_target_par(4,ntarget),j + + + ol_target=0.0_r8 + Nw=0.0_r8 + area_target_par=0.0_r8 + + do count=1,jall + i = weights_lgr_index_all(count) + ix = weights_eul_index_all1(count)!,1) + iy = weights_eul_index_all2(count)!,2) + ip = weights_eul_index_all3(count)!,3) + ! convert to 1D indexing of cubed-sphere + ii = (ip-1)*ncube*ncube+(iy-1)*ncube+ix! + wt = weights_all(count,1) + !determine terr_if + terr_if=0._r8 + if (terr(ii).GT.(1116.2-0.878*sgh_target(i))) terr_if=1. + ! (1): the lower left corner + ! (2): the lower right corner + ! (3): the upper right corner + ! (4): the upper left corner + !OL1 + if (lat_terr(ii) &!(ii)& + .GT.(target_corner_lat_deg(1,i)+target_center_lat(i))/2..and. & + lat_terr(ii) &!(ii)& + .LT.(target_corner_lat_deg(4,i)+target_center_lat(i))/2.) then + Nw(1,i)=Nw(1,i)+wt*terr_if + area_target_par(1,i)=area_target_par(1,i)+wt + endif + + !OL2 + if (lon_terr(ii) &!(ii)& + .GT.(target_corner_lon_deg(1,i)+target_center_lon(i))/2..and. & + lon_terr(ii) &!(ii)& + .LT.(target_corner_lon_deg(3,i)+target_center_lon(i))/2.) then + Nw(2,i)=Nw(2,i)+wt*terr_if + area_target_par(2,i)=area_target_par(2,i)+wt + end if + + + !OL3 + if (lon_terr(ii) &!(ii)& + .LT.target_center_lon(i).and. & + lat_terr(ii) &!(ii)& + .LT.target_center_lat(i).or. & + lon_terr(ii) &!(ii)& + .GT.target_center_lon(i).and. & + lat_terr(ii) &!(ii)& + .GT.target_center_lat(i)) then + Nw(3,i)=Nw(3,i)+wt*terr_if + area_target_par(3,i)=area_target_par(3,i)+wt + end if + + + !OL4 + if (lat_terr(ii) & !(ii)& + .GT.target_center_lat(i).and. & + lon_terr(ii) & !(ii)& + .LT.target_center_lon(i).or. & + lat_terr(ii) & !(ii)& + .LT.target_center_lat(i).and. & + lon_terr(ii) & !(ii)& + .GT.target_center_lon(i)) then + Nw(4,i)=Nw(4,i)+wt*terr_if + area_target_par(4,i)=area_target_par(4,i)+wt + end if + + !Nw(4,i)=Nw(4,i)+wt*terr_if + !area_target_par(4,i)=area_target_par(4,i)+wt + + + + do j=1,4 + ol_target(i,j)=Nw(j,i)/(area_target_par(j,i)+1e-14)!Nt(i)!/2.) + enddo + ol_target(i,:)=ol_target(i,:)*landfrac_target(i) + end do + where(abs(ol_target)<.001_r8.or.abs(ol_target).gt.1e+7) ol_target=0.0_r8 +end subroutine OLorig +!#endif +!===================== +!=================================================================== +!===================== +!#if 0 +subroutine OLgrid(terr,terrx,terry,wt,b,a,n,theta_in,hc,OLout) +!use physconst, only: rh2o,zvir,pi,rearth +!use abortutils +!Xie add +IMPLICIT NONE +integer,intent(in) :: n +real(r8),intent(in) :: hc,wt(n),terr(n),a,b,theta_in!a dy,b dx +real(r8),intent(in) :: terrx(n),terry(n) +real(r8),intent(out) :: OLout +real(r8) :: theta,theta1,theta2,rad,interval +real(r8) :: terr_count(n),terr_whole_count(n),cx(n),c1,c2 +!local +integer :: i +real(r8) :: j + terr_count=0.0_r8 + terr_whole_count=0.0_r8 + c1=0.0_r8 + c2=0.0_r8 + cx=0.0_r8 + !determine an acute theta in the triangle + !or minus 180 equilvalent acute angle + !then turn into radian + rad=4.0_r8*atan(1.0_r8)/180.0_r8 + interval=0.0_r8 + + !initialize + theta1=0.0_r8 + theta2=0.0_r8 + !set inside -360~360 + !this adds robustness of the scheme to different angle + theta1=MOD(theta_in,360._r8) + !set negative axis into 0~360 + if (theta1.ge.-360._r8.and.theta1.lt.0._r8) then + theta1=theta1+360._r8 + endif + !now we have only 0~360 angle + if (theta1.ge. 0._r8.and.theta1.le. 90._r8) then + theta=theta1*rad + theta2=theta1 + else if (theta1.gt. 90._r8.and.theta1.le. 180._r8) then + theta=(180._r8-theta1)*rad + theta2=180._r8-theta1 + else if (theta1.gt. 180._r8.and.theta1.le. 270._r8) then + theta=(theta1-180._r8)*rad + theta2=theta1-180._r8 + !we only use 0~180, so this makes similar to 1st and 2nd quadrant + else if (theta1.gt. 270._r8.and.theta1.le. 360._r8) then + theta=(360._r8-theta1)*rad + theta2=360._r8-theta1 + !we only use 0~180, so this makes similar to 1st and 2nd quadrant + endif + !we use theta2 to judge instead + !theta2=theta1 + !theta1=theta2 + !we restrict the angle in the first and second quadrant + !the third and fourth for OL are similar when theta is + !transformed by minus pi(180) + !two parallel lines are included + !xsin(theta)-ycos(theta)=c1 + !xsin(theta)-ycos(theta)=c2 + !xsin(theta)-ycos(theta)=cx,min(c1,c2) 0) .AND. (j < ncube_reconstruct)) THEN + beta = gp(j) + beta_next = gp(j+1) + ELSEIF (j == -1) THEN + beta = -piq - (gp(3) + piq) + beta_next = -piq - (gp(2) + piq) + ELSEIF (j == 0) THEN + beta = -piq - (gp(2) + piq) + beta_next = -piq + ELSEIF (j == ncube_reconstruct) THEN + beta = piq + beta_next = piq + (piq - gp(ncube_reconstruct-1)) + ELSEIF (j == ncube_reconstruct+1) THEN + beta = piq + (piq - gp(ncube_reconstruct-1)) + beta_next = piq + (piq - gp(ncube_reconstruct-2)) + ENDIF + + DO i = -1, ncube_reconstruct+1 + IF ((i > 0) .AND. (i < ncube_reconstruct)) THEN + alpha = gp(i) + alpha_next = gp(i+1) + ELSEIF (i == -1) THEN + alpha = -piq - (gp(3) + piq) + alpha_next = -piq - (gp(2) + piq) + ELSEIF (i == 0) THEN + alpha = -piq - (gp(2) + piq) + alpha_next = -piq + ELSEIF (i == ncube_reconstruct) THEN + alpha = piq + alpha_next = piq + (piq - gp(ncube_reconstruct-1)) + ELSEIF (i == ncube_reconstruct+1) THEN + alpha = piq + (piq - gp(ncube_reconstruct-1)) + alpha_next = piq + (piq - gp(ncube_reconstruct-2)) + ENDIF + abp_centroid(1,i,j) = & + I_10_ab(alpha_next,beta_next)-I_10_ab(alpha ,beta_next)+& + I_10_ab(alpha ,beta )-I_10_ab(alpha_next,beta ) +! - ASINH(COS(alpha_next) * TAN(beta_next)) & +! + ASINH(COS(alpha_next) * TAN(beta)) & +! + ASINH(COS(alpha) * TAN(beta_next)) & +! - ASINH(COS(alpha) * TAN(beta)) + + abp_centroid(2,i,j) = & + I_01_ab(alpha_next,beta_next)-I_01_ab(alpha ,beta_next)+& + I_01_ab(alpha ,beta )-I_01_ab(alpha_next,beta ) +! - ASINH(TAN(alpha_next) * COS(beta_next)) & +! + ASINH(TAN(alpha_next) * COS(beta)) & +! + ASINH(TAN(alpha) * COS(beta_next)) & +! - ASINH(TAN(alpha) * COS(beta)) + + !ADD PHL START + IF (order>2) THEN + ! TAN(alpha)^2 component + abp_centroid(3,i,j) =& + I_20_ab(alpha_next,beta_next)-I_20_ab(alpha ,beta_next)+& + I_20_ab(alpha ,beta )-I_20_ab(alpha_next,beta ) + + ! TAN(beta)^2 component + abp_centroid(4,i,j) = & + I_02_ab(alpha_next,beta_next)-I_02_ab(alpha ,beta_next)+& + I_02_ab(alpha ,beta )-I_02_ab(alpha_next,beta ) + + ! TAN(alpha) TAN(beta) component + abp_centroid(5,i,j) = & + I_11_ab(alpha_next,beta_next)-I_11_ab(alpha ,beta_next)+& + I_11_ab(alpha ,beta )-I_11_ab(alpha_next,beta ) + ENDIF + !ADD PHL END + ENDDO + ENDDO + +! +! PHL outcommented below +! + ! High order calculations +! IF (order > 2) THEN +! DO k = 1, nlon +! DO i = 1, int_nx(nlat,k)-1 +! IF ((int_itype(i,k) > 4) .AND. (int_np(1,i,k) == 1)) THEN +! abp_centroid(3, int_a(i,k), int_b(i,k)) = & +! abp_centroid(3, int_a(i,k), int_b(i,k)) + int_wt_2a(i,k) +! abp_centroid(4, int_a(i,k), int_b(i,k)) = & +! abp_centroid(4, int_a(i,k), int_b(i,k)) + int_wt_2b(i,k) +! abp_centroid(5, int_a(i,k), int_b(i,k)) = & +! abp_centroid(5, int_a(i,k), int_b(i,k)) + int_wt_2c(i,k) +! ENDIF +! ENDDO +! ENDDO +! ENDIF + + ! Normalize with element areas + DO j = -1, ncube_reconstruct+1 + IF ((j > 0) .AND. (j < ncube_reconstruct)) THEN + beta = gp(j) + beta_next = gp(j+1) + ELSEIF (j == -1) THEN + beta = -piq - (gp(3) + piq) + beta_next = -piq - (gp(2) + piq) + ELSEIF (j == 0) THEN + beta = -piq - (gp(2) + piq) + beta_next = -piq + ELSEIF (j == ncube_reconstruct) THEN + beta = piq + beta_next = piq + (piq - gp(ncube_reconstruct-1)) + ELSEIF (j == ncube_reconstruct+1) THEN + beta = piq + (piq - gp(ncube_reconstruct-1)) + beta_next = piq + (piq - gp(ncube_reconstruct-2)) + ENDIF + DO i = -1, ncube_reconstruct+1 + IF ((i > 0) .AND. (i < ncube_reconstruct)) THEN + alpha = gp(i) + alpha_next = gp(i+1) + ELSEIF (i == -1) THEN + alpha = -piq - (gp(3) + piq) + alpha_next = -piq - (gp(2) + piq) + ELSEIF (i == 0) THEN + alpha = -piq - (gp(2) + piq) + alpha_next = -piq + ELSEIF (i == ncube_reconstruct) THEN + alpha = piq + alpha_next = piq + (piq - gp(ncube_reconstruct-1)) + ELSEIF (i == ncube_reconstruct+1) THEN + alpha = piq + (piq - gp(ncube_reconstruct-1)) + alpha_next = piq + (piq - gp(ncube_reconstruct-2)) + ENDIF + + IF ((i > 0) .AND. (i < ncube_reconstruct) .AND. (j > 0) .AND. (j < ncube_reconstruct)) THEN + area = DAcube(i,j) + ELSE + area = EquiangularElementArea(alpha, alpha_next - alpha, & + beta, beta_next - beta) + ENDIF + + abp_centroid(1,i,j) = abp_centroid(1,i,j) / area + abp_centroid(2,i,j) = abp_centroid(2,i,j) / area + + IF (order > 2) THEN + IF ((i > 0) .AND. (i < ncube_reconstruct) .AND. (j > 0) .AND. (j < ncube_reconstruct)) THEN + abp_centroid(3,i,j) = abp_centroid(3,i,j) / area + abp_centroid(4,i,j) = abp_centroid(4,i,j) / area + abp_centroid(5,i,j) = abp_centroid(5,i,j) / area + ENDIF + ENDIF + ENDDO + ENDDO + + WRITE(*,*) '...Done computing ABP element centroids' + + END SUBROUTINE ComputeABPElementCentroids + +!------------------------------------------------------------------------------ +! FUNCTION EvaluateABPReconstruction +! +! Description: +! Evaluate the sub-grid scale reconstruction at the given point. +! +! Parameters: +! fcubehalo - Array of element values +! recons - Array of reconstruction coefficients +! a - Index of element in alpha direction (1 <= a <= ncube_reconstruct-1) +! b - Index of element in beta direction (1 <= b <= ncube_reconstruct-1) +! p - Panel index of element +! alpha - Alpha coordinate of evaluation point +! beta - Beta coordinate of evaluation point +! order - Order of the reconstruction +! value (OUT) - Result of function evaluation at given point +!------------------------------------------------------------------------------ + SUBROUTINE EvaluateABPReconstruction( & + fcubehalo, recons, a, b, p, alpha, beta, order, value) + IMPLICIT NONE + + ! Dummy variables + REAL (KIND=dbl_kind), DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6), & + INTENT(IN) :: fcubehalo + + REAL (KIND=dbl_kind), DIMENSION(:,:,:,:), INTENT(IN) :: recons + INTEGER (KIND=int_kind), INTENT(IN) :: a, b, p + REAL (KIND=dbl_kind), INTENT(IN) :: alpha, beta + INTEGER (KIND=int_kind), INTENT(IN) :: order + + REAL (KIND=dbl_kind), INTENT(OUT) :: value + + ! Evaluate constant order terms + value = fcubehalo(a,b,p) + + ! Evaluate linear order terms + IF (order > 1) THEN + value = value + recons(1,a,b,p) * (TAN(alpha) - abp_centroid(1,a,b)) + value = value + recons(2,a,b,p) * (TAN(beta) - abp_centroid(2,a,b)) + ENDIF + + ! Evaluate second order terms + IF (order > 2) THEN + value = value + recons(3,a,b,p) * & + (abp_centroid(1,a,b)**2 - abp_centroid(3,a,b)) + value = value + recons(4,a,b,p) * & + (abp_centroid(2,a,b)**2 - abp_centroid(4,a,b)) + value = value + recons(5,a,b,p) * & + (abp_centroid(1,a,b) * abp_centroid(2,a,b) - & + abp_centroid(5,a,b)) + + value = value + recons(3,a,b,p) * (TAN(alpha) - abp_centroid(1,a,b))**2 + value = value + recons(4,a,b,p) * (TAN(beta) - abp_centroid(2,a,b))**2 + value = value + recons(5,a,b,p) * (TAN(alpha) - abp_centroid(1,a,b)) & + * (TAN(beta) - abp_centroid(2,a,b)) + ENDIF + + END SUBROUTINE + +!------------------------------------------------------------------------------ +! SUBROUTINE ABPHaloMinMax +! +! Description: +! Calculate the minimum and maximum values of the cell-averaged function +! around the given element. +! +! Parameters: +! fcubehalo - Cell-averages for the cubed sphere +! a - Local element alpha index +! b - Local element beta index +! p - Local element panel index +! min_val (OUT) - Minimum value in the halo +! max_val (OUT) - Maximum value in the halo +! nomiddle - whether to not include the middle cell (index a,b) in the search. +! +! NOTE: Since this routine is not vectorized, it will likely be called MANY times. +! To speed things up, make sure to pass the first argument as the ENTIRE original +! array, not as a subset of it, since repeatedly cutting up that array and creating +! an array temporary (on some compilers) is VERY slow. +! ex: +! CALL APBHaloMinMax(zarg, a, ...) !YES +! CALL ABPHaloMinMax(zarg(-1:ncube_reconstruct+1,-1:ncube_reconstruct+1,:)) !NO -- slow +!------------------------------------------------------------------------------ + SUBROUTINE ABPHaloMinMax(fcubehalo, a, b, p, min_val, max_val, nomiddle) + IMPLICIT NONE + + REAL (KIND=dbl_kind), DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6), & + INTENT(IN) :: fcubehalo + + INTEGER (KIND=int_kind), INTENT(IN) :: a, b, p + REAL (KIND=dbl_kind), INTENT(OUT) :: min_val, max_val + LOGICAL, INTENT(IN) :: nomiddle + + ! Local variables + INTEGER (KIND=int_kind) :: i, j, il, jl, inew, jnew + REAL (KIND=dbl_kind) :: value + + min_val = fcubehalo(a,b,p) + max_val = fcubehalo(a,b,p) + value = fcubehalo(a,b,p) + + DO il = a-1,a+1 + DO jl = b-1,b+1 + + i = il + j = jl + + inew = i + jnew = j + + IF (nomiddle .AND. i==a .AND. j==b) CYCLE + + !Interior + IF ((i > 0) .AND. (i < ncube_reconstruct) .AND. (j > 0) .AND. (j < ncube_reconstruct)) THEN + value = fcubehalo(i,j,p) + + ELSE + + + !The next 4.0 regions are cases in which a,b themselves lie in the panel's halo, and the cell's "halo" (in this usage the 8.0 cells surrounding it) might wrap around into another part of the halo. This happens for (a,b) = {(1,:0),(ncube_reconstruct-1,:0),(1,ncube_reconstruct:),(ncube_reconstruct-1,ncube_reconstruct:)} and for the transposes thereof ({(:0,1), etc.}). In these cases (i,j) could lie in the "Corners" where nothing should lie. We correct this by moving i,j to its appropriate position on the "facing" halo, and then the remainder of the routine then moves it onto the correct face. + +101 FORMAT("ERROR cannot find (i,j) = (", I4, ", ", I4, ") for (a,b,p) = ", I4, ",", I4, ",", I4, ")") +102 FORMAT("i,j,p = ", 3I4, " moved to " 2I4, " (CASE ", I1, ")") + !NOTE: we need the general case to be able to properly handle (0,0), (ncube_reconstruct,0), etc. Note that we don't need to bother with (0,0), etc. when a, b lie in the interior, since both sides of the (0,0) cell are already accounted for by this routine. + !LOWER LEFT + IF (i < 1 .AND. j < 1) THEN + IF (a < 1) THEN !(a,b) centered on left halo, cross to lower halo + inew = 1-j + jnew = i + ELSE IF (b < 1) THEN !(a,b) centered on lower halo, cross to left halo + jnew = 1-i + inew = j + END IF +! WRITE(*,102) i, j, p, inew, jnew, 1 + !LOWER RIGHT + ELSE IF (i > ncube_reconstruct-1 .AND. j < 1) THEN + IF (a > ncube_reconstruct-1) THEN !(a,b) centered on right halo, cross to lower halo + inew = ncube_reconstruct-1+j + jnew = ncube_reconstruct-i + ELSE IF (b < 1) THEN !(a,b) centered on lower halo, cross to right halo + jnew = 1+(i-ncube_reconstruct) + inew = ncube_reconstruct-j + END IF +! WRITE(*,102) i, j, p, inew, jnew, 2 + !UPPER LEFT + ELSE IF (i < 1 .AND. j > ncube_reconstruct-1) THEN + IF (a < 1) THEN! (a,b) centered on left halo, cross to upper halo + inew = 1-(j-ncube_reconstruct) + jnew = ncube_reconstruct-i + ELSE IF (b > ncube_reconstruct-1) THEN !(a,b) centered on upper halo, cross to left halo + inew = ncube_reconstruct-j + jnew = ncube_reconstruct-1-i + END IF +! WRITE(*,102) i, j, p, inew, jnew, 3 + !UPPER RIGHT + ELSE IF (i > ncube_reconstruct-1 .AND. j > ncube_reconstruct-1) THEN + IF (a > ncube_reconstruct-1) THEN !(a,b) centered on right halo, cross to upper halo + inew = ncube_reconstruct-1-(ncube_reconstruct-j) + jnew = i + ELSE IF (b > ncube_reconstruct-1) THEN !(a,b) centered on upper halo, cross to right halo + inew = j + jnew = ncube_reconstruct-1-(ncube_reconstruct-i) + END IF +! WRITE(*,102) i, j, p, inew, jnew, 4 + END IF + + i = inew + j = jnew + + + !Lower halo ("halo" meaning the panel's halo, not the nine-cell halo + IF ((i < 1) .AND. (j > 0) .AND. (j < ncube_reconstruct)) THEN + IF (p == 1) THEN + value = fcubehalo(ncube_reconstruct-1+i,j,4) + ELSEIF (p == 2) THEN + value = fcubehalo(ncube_reconstruct-1+i,j,1) + ELSEIF (p == 3) THEN + value = fcubehalo(ncube_reconstruct-1+i,j,2) + ELSEIF (p == 4) THEN + value = fcubehalo(ncube_reconstruct-1+i,j,3) + ELSEIF (p == 5) THEN + value = fcubehalo(j,1-i,4) + ELSEIF (p == 6) THEN + value = fcubehalo(ncube_reconstruct-j,ncube_reconstruct-1+i,4) + ENDIF + + !Upper halo + ELSEIF ((i > ncube_reconstruct-1) .AND. (j > 0) .AND. (j < ncube_reconstruct)) THEN + IF (p == 1) THEN + value = fcubehalo(i-ncube_reconstruct+1,j,2) + ELSEIF (p == 2) THEN + value = fcubehalo(i-ncube_reconstruct+1,j,3) + ELSEIF (p == 3) THEN + value = fcubehalo(i-ncube_reconstruct+1,j,4) + ELSEIF (p == 4) THEN + value = fcubehalo(i-ncube_reconstruct+1,j,1) + ELSEIF (p == 5) THEN + value = fcubehalo(ncube_reconstruct-j,i-ncube_reconstruct+1,2) + ELSEIF (p == 6) THEN + value = fcubehalo(j,2*ncube_reconstruct-i-1,2) + ENDIF + + !Left halo + ELSEIF ((j < 1) .AND. (i > 0) .AND. (i < ncube_reconstruct)) THEN + IF (p == 1) THEN + value = fcubehalo(i,ncube_reconstruct-1+j,5) + ELSEIF (p == 2) THEN + value = fcubehalo(ncube_reconstruct-1+j,ncube_reconstruct-i,5) + ELSEIF (p == 3) THEN + value = fcubehalo(ncube_reconstruct-i,1-j,5) + ELSEIF (p == 4) THEN + value = fcubehalo(1-j,i,5) + ELSEIF (p == 5) THEN + value = fcubehalo(ncube_reconstruct-i,1-j,3) + ELSEIF (p == 6) THEN + value = fcubehalo(i,ncube_reconstruct-1+j,1) + ENDIF + + !Right halo + ELSEIF ((j > ncube_reconstruct-1) .AND. (i > 0) .AND. (i < ncube_reconstruct)) THEN + IF (p == 1) THEN + value = fcubehalo(i,j-ncube_reconstruct+1,6) + ELSEIF (p == 2) THEN + value = fcubehalo(2*ncube_reconstruct-j-1,i,6) + ELSEIF (p == 3) THEN + value = fcubehalo(ncube_reconstruct-i, 2*ncube_reconstruct-j-1,6) + ELSEIF (p == 4) THEN + value = fcubehalo(j-ncube_reconstruct+1,ncube_reconstruct-i,6) + ELSEIF (p == 5) THEN + value = fcubehalo(i,j-ncube_reconstruct+1,1) + ELSEIF (p == 6) THEN + value = fcubehalo(ncube_reconstruct-i, 2*ncube_reconstruct-j-1,3) + ENDIF + + ENDIF + + END IF + min_val = MIN(min_val, value) + max_val = MAX(max_val, value) + ENDDO + ENDDO + END SUBROUTINE + +!------------------------------------------------------------------------------ +! SUBROUTINE MonotonizeABPGradient +! +! Description: +! Apply a monotonic filter to the calculated ABP gradient. +! +! Parameters: +! fcubehalo - Scalar field on the cubed sphere to use in reconstruction +! order - Order of the reconstruction +! recons (INOUT) - Array of reconstructed coefficients +! selective - whether to apply a simple form of selective limiting, + !which assumes that if a point is larger/smaller than ALL of its + !surrounding points, that the extremum is physical, and that + !filtering should not be applied to it. +! +! Remarks: +! This monotonizing scheme is based on the monotone scheme for unstructured +! grids of Barth and Jespersen (1989). +!------------------------------------------------------------------------------ + SUBROUTINE MonotonizeABPGradient(fcubehalo, order, recons, selective) + +! USE selective_limiting + + IMPLICIT NONE + + REAL (KIND=dbl_kind), DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6), & + INTENT(IN) :: fcubehalo + + INTEGER (KIND=int_kind), INTENT(IN) :: order + + LOGICAL, INTENT(IN) :: selective + + REAL (KIND=dbl_kind), DIMENSION(:,:,:,:), INTENT(INOUT) :: recons + + ! Local variables + INTEGER (KIND=int_kind) :: i, j, k, m, n, skip + + REAL (KIND=dbl_kind) :: local_min, local_max, value, phi, min_phi + REAL (KIND=dbl_kind) :: disc, mx, my, lam, gamma_min, gamma_max + REAL (KIND=dbl_kind), DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6) :: & + gamma + + ! The first-order piecewise constant scheme is monotone by construction + IF (order == 1) THEN + RETURN + ENDIF + +! +! xxxxx +! +! IF (selective) THEN +! CALL smoothness2D(fcubehalo, gamma, 2) +! WRITE(*,*) 'gamma range: max ', MAXVAL(gamma), " min ", MINVAL(gamma) +! DO i=1,ncube_reconstruct-1 +! WRITE(*,*) gamma(i, i, 3) +! ENDDO +! skip = 0 +! END IF + + + ! Apply monotone limiting + DO k = 1, 6 + DO j = 1, ncube_reconstruct-1 + DO i = 1, ncube_reconstruct-1 + + + IF (selective) THEN + + CALL ABPHaloMinMax(gamma, i, j, k, gamma_min, gamma_max, .FALSE.) + + IF (gamma_max/(gamma_min + tiny) < lammax) THEN + skip = skip + 1 + CYCLE + END IF + + END IF + + CALL ABPHaloMinMax(fcubehalo, i, j, k, local_min, local_max,.FALSE.) + + + ! Initialize the limiter + min_phi = one + + ! For the second-order calculation, the minima and maxima will occur + ! at the corner points of the element + DO m = i, i+1 + DO n = j, j+1 + + ! Evaluate the function at each corner point + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, gp(m), gp(n), order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), local_min, local_max, min_phi) + ENDDO + ENDDO + + ! For the third order method, the minima and maxima may occur along + ! the line segments given by du/dx = 0 and du/dy = 0. Also check + ! for the presence of a maxima / minima of the quadratic within + ! the domain. + IF (order == 3) THEN + disc = recons(5,i,j,k)**2 - 4.0 * recons(4,i,j,k) * recons(3,i,j,k) + + ! Check if the quadratic is minimized within the element + IF (ABS(disc) > tiny) THEN + mx = - recons(5,i,j,k) * recons(2,i,j,k) & + + 2.0 * recons(4,i,j,k) * recons(1,i,j,k) + my = - recons(5,i,j,k) * recons(1,i,j,k) & + + 2.0 * recons(3,i,j,k) * recons(2,i,j,k) + + mx = mx / disc + abp_centroid(1,i,j) + my = my / disc + abp_centroid(2,i,j) + + IF ((mx - TAN(gp(i)) > -tiny) .AND. & + (mx - TAN(gp(i+1)) < tiny) .AND. & + (my - TAN(gp(j)) > -tiny) .AND. & + (my - TAN(gp(j+1)) < tiny) & + ) THEN + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, ATAN(mx), ATAN(my), & + order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), & + local_min, local_max, min_phi) + ENDIF + ENDIF + + ! Check all potential minimizer points along element boundaries + IF (ABS(recons(5,i,j,k)) > tiny) THEN + + ! Left/right edge, intercept with du/dx = 0 + DO m = i, i+1 + my = - recons(1,i,j,k) - 2.0 * recons(3,i,j,k) * & + (TAN(gp(m)) - abp_centroid(1,i,j)) + + my = my / recons(5,i,j,k) + abp_centroid(2,i,j) + + IF ((my < TAN(gp(j))) .OR. (my > TAN(gp(j+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, gp(m), ATAN(my), & + order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), & + local_min, local_max, min_phi) + ENDDO + + ! Top/bottom edge, intercept with du/dy = 0 + DO n = j, j+1 + mx = - recons(2,i,j,k) - 2.0 * recons(4,i,j,k) * & + (TAN(gp(n)) - abp_centroid(2,i,j)) + + mx = mx / recons(5,i,j,k) + abp_centroid(1,i,j) + + IF ((mx < TAN(gp(i))) .OR. (mx > TAN(gp(i+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, ATAN(mx), gp(n), & + order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), & + local_min, local_max, min_phi) + ENDDO + ENDIF + + ! Top/bottom edge, intercept with du/dx = 0 + IF (ABS(recons(3,i,j,k)) > tiny) THEN + DO n = j, j+1 + mx = - recons(1,i,j,k) - recons(5,i,j,k) * & + (TAN(gp(n)) - abp_centroid(2,i,j)) + + mx = mx / (2.0 * recons(3,i,j,k)) + abp_centroid(1,i,j) + + IF ((mx < TAN(gp(i))) .OR. (mx > TAN(gp(i+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, ATAN(mx), gp(n), & + order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), & + local_min, local_max, min_phi) + ENDDO + ENDIF + + ! Left/right edge, intercept with du/dy = 0 + IF (ABS(recons(4,i,j,k)) > tiny) THEN + DO m = i, i+1 + my = - recons(2,i,j,k) - recons(5,i,j,k) * & + (TAN(gp(m)) - abp_centroid(1,i,j)) + + my = my / (2.0 * recons(4,i,j,k)) + abp_centroid(2,i,j) + + IF ((my < TAN(gp(j))) .OR. (my > TAN(gp(j+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, gp(m), ATAN(my), & + order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), & + local_min, local_max, min_phi) + ENDDO + ENDIF + ENDIF + + IF ((min_phi < -tiny) .OR. (min_phi > one + tiny)) THEN + WRITE (*,*) 'Fatal Error: In MonotonizeABPGradient' + WRITE (*,*) 'Slope limiter out of range: ', min_phi + STOP + ENDIF + + ! Apply monotone limiter to all reconstruction coefficients + recons(1,i,j,k) = min_phi * recons(1,i,j,k) + recons(2,i,j,k) = min_phi * recons(2,i,j,k) + + IF (order > 2) THEN + recons(3,i,j,k) = min_phi * recons(3,i,j,k) + recons(4,i,j,k) = min_phi * recons(4,i,j,k) + recons(5,i,j,k) = min_phi * recons(5,i,j,k) + ENDIF + ENDDO + ENDDO + ENDDO + + IF (selective) WRITE(*,*) 'skipped ', skip, ' points out of ', 6*(ncube_reconstruct-1)**2 + + END SUBROUTINE + +!------------------------------------------------------------------------------ +! SUBROUTINE PosDefABPGradient +! +! Description: +! Scale the reconstructions so they are positive definite +! +! Parameters: +! fcubehalo - Scalar field on the cubed sphere to use in reconstruction +! order - Order of the reconstruction +! recons (INOUT) - Array of reconstructed coefficients +! +! Remarks: +! This monotonizing scheme is based on the monotone scheme for unstructured +! grids of Barth and Jespersen (1989), but simpler. This simply finds the +! minimum and then scales the reconstruction so that it is 0. +!------------------------------------------------------------------------------ + SUBROUTINE PosDefABPGradient(fcubehalo, order, recons) + + IMPLICIT NONE + + REAL (KIND=dbl_kind), DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6), & + INTENT(IN) :: fcubehalo + + INTEGER (KIND=int_kind), INTENT(IN) :: order + + REAL (KIND=dbl_kind), DIMENSION(:,:,:,:), INTENT(INOUT) :: recons + + ! Local variables + INTEGER (KIND=int_kind) :: i, j, k, m, n + + REAL (KIND=dbl_kind) :: local_min, local_max, value, phi, min_phi + REAL (KIND=dbl_kind) :: disc, mx, my, lam, gamma_min, gamma_max + REAL (KIND=dbl_kind), DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6) :: & + gamma + + ! The first-order piecewise constant scheme is monotone by construction + IF (order == 1) THEN + RETURN + ENDIF + + + ! Apply monotone limiting + DO k = 1, 6 + DO j = 1, ncube_reconstruct-1 + DO i = 1, ncube_reconstruct-1 + + !If the average value in the cell is 0.0, then we should skip + !all of the scaling and just set the reconstruction to 0.0 +! IF (ABS(fcubehalo(i,j,k)) < tiny) THEN +! recons(:,i,j,k) = 0.0 +! CYCLE +! END IF + CALL ABPHaloMinMax(fcubehalo, i, j, k, local_min, local_max,.FALSE.) + + + !This allowance for miniscule negative values appearing around the cell being + !filtered/limited. Before this, negative values would be caught in adjust_limiter + !and would stop the model. Doing this only causes minor negative values; no blowing + !up is observed. The rationale is the same as for the monotone filter, which does + !allow miniscule negative values due to roundoff error --- of the order E-10 --- + !in flux-form methods (and E-17 in the s-L method, indicating that roundoff error + !is more severe in the flux-form method, as we expect since we are often subtracting + !2.0 values which are very close together. + local_min = MIN(0.0,local_min) + local_max = bignum !prevents scaling upward; for positive + !definite limiting we don't care about the upper bound + + ! Initialize the limiter + min_phi = one + + ! For the second-order calculation, the minima and maxima will occur + ! at the corner points of the element + DO m = i, i+1 + DO n = j, j+1 + + ! Evaluate the function at each corner point + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, gp(m), gp(n), order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), local_min, local_max, min_phi) + ENDDO + ENDDO + + ! For the third order method, the minima and maxima may occur along + ! the line segments given by du/dx = 0 and du/dy = 0. Also check + ! for the presence of a maxima / minima of the quadratic within + ! the domain. + IF (order == 3) THEN + disc = recons(5,i,j,k)**2 - 4.0 * recons(4,i,j,k) * recons(3,i,j,k) + + ! Check if the quadratic is minimized within the element + IF (ABS(disc) > tiny) THEN + mx = - recons(5,i,j,k) * recons(2,i,j,k) & + + 2.0 * recons(4,i,j,k) * recons(1,i,j,k) + my = - recons(5,i,j,k) * recons(1,i,j,k) & + + 2.0 * recons(3,i,j,k) * recons(2,i,j,k) + + mx = mx / disc + abp_centroid(1,i,j) + my = my / disc + abp_centroid(2,i,j) + + IF ((mx - TAN(gp(i)) > -tiny) .AND. & + (mx - TAN(gp(i+1)) < tiny) .AND. & + (my - TAN(gp(j)) > -tiny) .AND. & + (my - TAN(gp(j+1)) < tiny) & + ) THEN + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, ATAN(mx), ATAN(my), & + order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), & + local_min, local_max, min_phi) + ENDIF + ENDIF + + ! Check all potential minimizer points along element boundaries + IF (ABS(recons(5,i,j,k)) > tiny) THEN + + ! Left/right edge, intercept with du/dx = 0 + DO m = i, i+1 + my = - recons(1,i,j,k) - 2.0 * recons(3,i,j,k) * & + (TAN(gp(m)) - abp_centroid(1,i,j)) + + my = my / recons(5,i,j,k) + abp_centroid(2,i,j) + + IF ((my < TAN(gp(j))) .OR. (my > TAN(gp(j+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, gp(m), ATAN(my), & + order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), & + local_min, local_max, min_phi) + ENDDO + + ! Top/bottom edge, intercept with du/dy = 0 + DO n = j, j+1 + mx = - recons(2,i,j,k) - 2.0 * recons(4,i,j,k) * & + (TAN(gp(n)) - abp_centroid(2,i,j)) + + mx = mx / recons(5,i,j,k) + abp_centroid(1,i,j) + + IF ((mx < TAN(gp(i))) .OR. (mx > TAN(gp(i+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, ATAN(mx), gp(n), & + order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), & + local_min, local_max, min_phi) + ENDDO + ENDIF + + ! Top/bottom edge, intercept with du/dx = 0 + IF (ABS(recons(3,i,j,k)) > tiny) THEN + DO n = j, j+1 + mx = - recons(1,i,j,k) - recons(5,i,j,k) * & + (TAN(gp(n)) - abp_centroid(2,i,j)) + + mx = mx / (2.0 * recons(3,i,j,k)) + abp_centroid(1,i,j) + + IF ((mx < TAN(gp(i))) .OR. (mx > TAN(gp(i+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, ATAN(mx), gp(n), & + order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), & + local_min, local_max, min_phi) + ENDDO + ENDIF + + ! Left/right edge, intercept with du/dy = 0 + IF (ABS(recons(4,i,j,k)) > tiny) THEN + DO m = i, i+1 + my = - recons(2,i,j,k) - recons(5,i,j,k) * & + (TAN(gp(m)) - abp_centroid(1,i,j)) + + my = my / (2.0 * recons(4,i,j,k)) + abp_centroid(2,i,j) + + IF ((my < TAN(gp(j))) .OR. (my > TAN(gp(j+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, gp(m), ATAN(my), & + order, value) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), & + local_min, local_max, min_phi) + ENDDO + ENDIF + ENDIF + + IF ((min_phi < -tiny) .OR. (min_phi > one + tiny)) THEN + WRITE (*,*) 'Fatal Error: In MonotonizeABPGradient' + WRITE (*,*) 'Slope limiter out of range: ', min_phi + STOP + ENDIF + + ! Apply monotone limiter to all reconstruction coefficients + recons(1,i,j,k) = min_phi * recons(1,i,j,k) + recons(2,i,j,k) = min_phi * recons(2,i,j,k) + + IF (order > 2) THEN + recons(3,i,j,k) = min_phi * recons(3,i,j,k) + recons(4,i,j,k) = min_phi * recons(4,i,j,k) + recons(5,i,j,k) = min_phi * recons(5,i,j,k) + ENDIF + + ENDDO + ENDDO + ENDDO + + + END SUBROUTINE PosDefABPGradient + +!------------------------------------------------------------------------------ +! SUBROUTINE MonotonizeABPGradient_New +! +! Description: +! Apply a monotonic filter to the calculated ABP gradient. +! +! Parameters: +! fcubehalo - Scalar field on the cubed sphere to use in reconstruction +! order - Order of the reconstruction +! recons (INOUT) - Array of reconstructed coefficients +! +! Remarks: +! This monotonizing scheme is similar to the one in MonotonizeABPGradient, +! except the second order derivatives are limited after the first order +! derivatives. +!------------------------------------------------------------------------------ + SUBROUTINE MonotonizeABPGradient_New(fcubehalo, order, recons) + + IMPLICIT NONE + + REAL (KIND=dbl_kind), DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6), & + INTENT(IN) :: fcubehalo + + INTEGER (KIND=int_kind), INTENT(IN) :: order + + REAL (KIND=dbl_kind), DIMENSION(:,:,:,:), INTENT(INOUT) :: recons + + ! Local variables + INTEGER (KIND=int_kind) :: i, j, k, m, n + + REAL (KIND=dbl_kind) :: local_min, local_max, value, phi, min_phi, linval + REAL (KIND=dbl_kind) :: disc, mx, my + + ! The first-order piecewise constant scheme is monotone by construction + IF (order == 1) THEN + RETURN + ENDIF + + ! Apply monotone limiting + DO k = 1, 6 + DO j = 1, ncube_reconstruct-1 + DO i = 1, ncube_reconstruct-1 + CALL ABPHaloMinMax(fcubehalo, i, j, k, local_min, local_max, .FALSE.) + + ! Initialize the limiter + min_phi = one + + ! For the second-order calculation, the minima and maxima will occur + ! at the corner points of the element + DO m = i, i+1 + DO n = j, j+1 + + ! Evaluate the function at each corner point, only taking into + ! account the linear component of the reconstruction. + value = & + fcubehalo(i,j,k) & + + recons(1,i,j,k) * (TAN(gp(m)) - abp_centroid(1,i,j)) & + + recons(2,i,j,k) * (TAN(gp(n)) - abp_centroid(2,i,j)) + + CALL AdjustLimiter( & + value, fcubehalo(i,j,k), local_min, local_max, min_phi) + ENDDO + ENDDO + + ! Apply monotone limiter to all reconstruction coefficients + IF ((min_phi < -tiny) .OR. (min_phi > one + tiny)) THEN + WRITE (*,*) 'Fatal Error: In MonotonizeABPGradient' + WRITE (*,*) 'Slope limiter out of range: ', min_phi + STOP + ENDIF + + recons(1,i,j,k) = min_phi * recons(1,i,j,k) + recons(2,i,j,k) = min_phi * recons(2,i,j,k) + + ! For the third order method, the minima and maxima may occur along + ! the line segments given by du/dx = 0 and du/dy = 0. Also check + ! for the presence of a maxima / minima of the quadratic within + ! the domain. + IF (order == 3) THEN + ! Reset the limiter + min_phi = one + + ! Calculate discriminant, which we use to determine the absolute + ! minima/maxima of the paraboloid + disc = recons(5,i,j,k)**2 - 4.0 * recons(4,i,j,k) * recons(3,i,j,k) + + ! Check if the quadratic is minimized within the element + IF (ABS(disc) > tiny) THEN + mx = - recons(5,i,j,k) * recons(2,i,j,k) & + + 2.0 * recons(4,i,j,k) * recons(1,i,j,k) + my = - recons(5,i,j,k) * recons(1,i,j,k) & + + 2.0 * recons(3,i,j,k) * recons(2,i,j,k) + + mx = mx / disc + abp_centroid(1,i,j) + my = my / disc + abp_centroid(2,i,j) + + IF ((mx - TAN(gp(i)) > -tiny) .AND. & + (mx - TAN(gp(i+1)) < tiny) .AND. & + (my - TAN(gp(j)) > -tiny) .AND. & + (my - TAN(gp(j+1)) < tiny) & + ) THEN + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, ATAN(mx), ATAN(my), & + order, value) + + linval = & + fcubehalo(i,j,k) & + + recons(1,i,j,k) * (mx - abp_centroid(1,i,j)) & + + recons(2,i,j,k) * (my - abp_centroid(2,i,j)) + + IF (linval < local_min) THEN + linval = local_min + ENDIF + IF (linval > local_max) THEN + linval = local_max + ENDIF + + CALL AdjustLimiter( & + value, linval, local_min, local_max, min_phi) + ENDIF + ENDIF + + ! Check all potential minimizer points along element boundaries + IF (ABS(recons(5,i,j,k)) > tiny) THEN + + ! Left/right edge, intercept with du/dx = 0 + DO m = i, i+1 + my = - recons(1,i,j,k) - 2.0 * recons(3,i,j,k) * & + (TAN(gp(m)) - abp_centroid(1,i,j)) + + my = my / recons(5,i,j,k) + abp_centroid(2,i,j) + + IF ((my < TAN(gp(j))) .OR. (my > TAN(gp(j+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, gp(m), ATAN(my), & + order, value) + + linval = & + fcubehalo(i,j,k) & + + recons(1,i,j,k) * (TAN(gp(m)) - abp_centroid(1,i,j)) & + + recons(2,i,j,k) * (my - abp_centroid(2,i,j)) + + IF (linval < local_min) THEN + linval = local_min + ENDIF + IF (linval > local_max) THEN + linval = local_max + ENDIF + + CALL AdjustLimiter( & + value, linval, local_min, local_max, min_phi) + ENDDO + + ! Top/bottom edge, intercept with du/dy = 0 + DO n = j, j+1 + mx = - recons(2,i,j,k) - 2.0 * recons(4,i,j,k) * & + (TAN(gp(n)) - abp_centroid(2,i,j)) + + mx = mx / recons(5,i,j,k) + abp_centroid(1,i,j) + + IF ((mx < TAN(gp(i))) .OR. (mx > TAN(gp(i+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, ATAN(mx), gp(n), & + order, value) + + linval = & + fcubehalo(i,j,k) & + + recons(1,i,j,k) * (mx - abp_centroid(1,i,j)) & + + recons(2,i,j,k) * (TAN(gp(n)) - abp_centroid(2,i,j)) + + IF (linval < local_min) THEN + linval = local_min + ENDIF + IF (linval > local_max) THEN + linval = local_max + ENDIF + + CALL AdjustLimiter( & + value, linval, local_min, local_max, min_phi) + ENDDO + ENDIF + + ! Top/bottom edge, intercept with du/dx = 0 + IF (ABS(recons(3,i,j,k)) > tiny) THEN + DO n = j, j+1 + mx = - recons(1,i,j,k) - recons(5,i,j,k) * & + (TAN(gp(n)) - abp_centroid(2,i,j)) + + mx = mx / (2.0 * recons(3,i,j,k)) + abp_centroid(1,i,j) + + IF ((mx < TAN(gp(i))) .OR. (mx > TAN(gp(i+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, ATAN(mx), gp(n), & + order, value) + + linval = & + fcubehalo(i,j,k) & + + recons(1,i,j,k) * (mx - abp_centroid(1,i,j)) & + + recons(2,i,j,k) * (TAN(gp(n)) - abp_centroid(2,i,j)) + + IF (linval < local_min) THEN + linval = local_min + ENDIF + IF (linval > local_max) THEN + linval = local_max + ENDIF + + CALL AdjustLimiter( & + value, linval, local_min, local_max, min_phi) + ENDDO + ENDIF + + ! Left/right edge, intercept with du/dy = 0 + IF (ABS(recons(4,i,j,k)) > tiny) THEN + DO m = i, i+1 + my = - recons(2,i,j,k) - recons(5,i,j,k) * & + (TAN(gp(m)) - abp_centroid(1,i,j)) + + my = my / (2.0 * recons(4,i,j,k)) + abp_centroid(2,i,j) + + IF ((my < TAN(gp(j))) .OR. (my > TAN(gp(j+1)))) THEN + CYCLE + ENDIF + + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, gp(m), ATAN(my), & + order, value) + + linval = & + fcubehalo(i,j,k) & + + recons(1,i,j,k) * (TAN(gp(m)) - abp_centroid(1,i,j)) & + + recons(2,i,j,k) * (my - abp_centroid(2,i,j)) + + IF (linval < local_min) THEN + linval = local_min + ENDIF + IF (linval > local_max) THEN + linval = local_max + ENDIF + + CALL AdjustLimiter( & + value, linval, local_min, local_max, min_phi) + ENDDO + ENDIF + + ! For the second-order calculation, the minima and maxima will occur + ! at the corner points of the element + DO m = i, i+1 + DO n = j, j+1 + + ! Evaluate the function at each corner point + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, gp(m), gp(n), & + order, value) + + linval = & + fcubehalo(i,j,k) & + + recons(1,i,j,k) * (TAN(gp(m)) - abp_centroid(1,i,j)) & + + recons(2,i,j,k) * (TAN(gp(n)) - abp_centroid(2,i,j)) + + IF (linval < local_min) THEN + linval = local_min + ENDIF + IF (linval > local_max) THEN + linval = local_max + ENDIF + + CALL AdjustLimiter( & + value, linval, local_min, local_max, min_phi) + ENDDO + ENDDO + + IF ((min_phi < -tiny) .OR. (min_phi > one + tiny)) THEN + WRITE (*,*) 'Fatal Error: In MonotonizeABPGradient' + WRITE (*,*) 'Slope limiter out of range: ', min_phi + STOP + ENDIF + + WRITE (*,*) '2: ', min_phi + + recons(1,i,j,k) = min_phi * recons(1,i,j,k) + recons(2,i,j,k) = min_phi * recons(2,i,j,k) + recons(3,i,j,k) = min_phi * recons(3,i,j,k) + recons(4,i,j,k) = min_phi * recons(4,i,j,k) + recons(5,i,j,k) = min_phi * recons(5,i,j,k) + ENDIF + ENDDO + ENDDO + ENDDO + + END SUBROUTINE + +!------------------------------------------------------------------------------ +! SUBROUTINE ReconstructABPGradient_NEL +! +! Description: +! Construct a non-equidistant linear reconstruction of the gradient +! within each element on an ABP grid. +! +! Parameters: +! fcubehalo - Scalar field on the ABP grid to use in reconstruction +! recons (OUT) - Array of reconstructed coefficients for total elements +! order - Order of the scheme (2 or 3) +!------------------------------------------------------------------------------ + SUBROUTINE ReconstructABPGradient_NEL(fcubehalo, recons, order) + +! USE CubedSphereTrans +! USE InterpolateCSLL_Utils + + IMPLICIT NONE + + REAL (KIND=dbl_kind), & + DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6), INTENT(IN) :: fcubehalo + + REAL (KIND=dbl_kind), DIMENSION(:,:,:,:), INTENT(OUT) :: recons + + INTEGER (KIND=int_kind), INTENT(IN) :: order + + ! Local variables + INTEGER (KIND=int_kind) :: i, j, p + + REAL (KIND=dbl_kind) :: alpha1, alpha2, beta1, beta2 + REAL (KIND=dbl_kind) :: dx_left, dx_right, top_value, bot_value + + DO p = 1, 6 + DO j = 1, ncube_reconstruct-1 + DO i = 1, ncube_reconstruct-1 + dx_left = abp_centroid(1,i-1,j) - abp_centroid(1,i,j) + dx_right = abp_centroid(1,i+1,j) - abp_centroid(1,i,j) + + recons(1,i,j,p) = & + (+ fcubehalo(i-1,j,p) * dx_right**2 & + - fcubehalo(i+1,j,p) * dx_left**2 & + - fcubehalo(i,j,p) * (dx_right**2 - dx_left**2)) / & + (dx_right * dx_left * (dx_right - dx_left)) + + dx_left = abp_centroid(2,i,j-1) - abp_centroid(2,i,j) + dx_right = abp_centroid(2,i,j+1) - abp_centroid(2,i,j) + + recons(2,i,j,p) = & + (+ fcubehalo(i,j-1,p) * dx_right**2 & + - fcubehalo(i,j+1,p) * dx_left**2 & + - fcubehalo(i,j,p) * (dx_right**2 - dx_left**2)) / & + (dx_right * dx_left * (dx_right - dx_left)) + + IF (order > 2) THEN + dx_left = abp_centroid(1,i-1,j) - abp_centroid(1,i,j) + dx_right = abp_centroid(1,i+1,j) - abp_centroid(1,i,j) + + recons(3,i,j,p) = & + (+ fcubehalo(i-1,j,p) * dx_right & + - fcubehalo(i+1,j,p) * dx_left & + - fcubehalo(i,j,p) * (dx_right - dx_left)) / & + (dx_right * dx_left * (dx_left - dx_right)) + + dx_left = abp_centroid(2,i,j-1) - abp_centroid(2,i,j) + dx_right = abp_centroid(2,i,j+1) - abp_centroid(2,i,j) + + recons(4,i,j,p) = & + (+ fcubehalo(i,j-1,p) * dx_right & + - fcubehalo(i,j+1,p) * dx_left & + - fcubehalo(i,j,p) * (dx_right - dx_left)) / & + (dx_right * dx_left * (dx_left - dx_right)) + ENDIF + ENDDO + ENDDO + + IF (order > 2) THEN + DO j = 1, ncube_reconstruct-1 + DO i = 1, ncube_reconstruct-1 + dx_left = abp_centroid(1,i-1,j+1) - abp_centroid(1,i,j+1) + dx_right = abp_centroid(1,i+1,j+1) - abp_centroid(1,i,j+1) + + top_value = & + (+ fcubehalo(i-1,j+1,p) * dx_right**2 & + - fcubehalo(i+1,j+1,p) * dx_left**2 & + - fcubehalo(i,j+1,p) * (dx_right**2 - dx_left**2)) / & + (dx_right * dx_left * (dx_right - dx_left)) + + dx_left = abp_centroid(1,i-1,j-1) - abp_centroid(1,i,j-1) + dx_right = abp_centroid(1,i+1,j-1) - abp_centroid(1,i,j-1) + + bot_value = & + (+ fcubehalo(i-1,j-1,p) * dx_right**2 & + - fcubehalo(i+1,j-1,p) * dx_left**2 & + - fcubehalo(i,j-1,p) * (dx_right**2 - dx_left**2)) / & + (dx_right * dx_left * (dx_right - dx_left)) + + dx_left = abp_centroid(2,i,j-1) - abp_centroid(2,i,j) + dx_right = abp_centroid(2,i,j+1) - abp_centroid(2,i,j) + + recons(5,i,j,p) = & + (+ bot_value * dx_right**2 & + - top_value * dx_left**2 & + - recons(1,i,j,p) * (dx_right**2 - dx_left**2)) / & + (dx_right * dx_left * (dx_right - dx_left)) + + ENDDO + ENDDO + ENDIF + ENDDO + + END SUBROUTINE + +!------------------------------------------------------------------------------ +! SUBROUTINE ReconstructABPGradient_NEP +! +! Description: +! Construct a non-equidistant parabolic reconstruction of the gradient +! within each element on an ABP grid. +! +! Parameters: +! fcubehalo - Scalar field on the ABP grid to use in reconstruction +! recons (OUT) - Array of reconstructed coefficients for total elements +! order - Order of the scheme (2 or 3) +!------------------------------------------------------------------------------ + SUBROUTINE ReconstructABPGradient_NEP(fcubehalo, recons, order) + + +! USE CubedSphereTrans +! USE InterpolateCSLL_Utils + + IMPLICIT NONE + + REAL (KIND=dbl_kind), & + DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6), INTENT(IN) :: fcubehalo + + REAL (KIND=dbl_kind), DIMENSION(:,:,:,:), INTENT(OUT) :: recons + + INTEGER (KIND=int_kind), INTENT(IN) :: order + + ! Local variables + INTEGER (KIND=int_kind) :: i, j, p + + REAL (KIND=dbl_kind) :: x1, x2, x4, x5, y1, y2, y3, y4, y5 + + REAL (KIND=dbl_kind), DIMENSION(5) :: t, pa, denom + + DO p = 1, 6 + DO j = 1, ncube_reconstruct-1 + DO i = 1, ncube_reconstruct-1 + ! X-direction reconstruction + x1 = abp_centroid(1,i-2,j) - abp_centroid(1,i,j) + x2 = abp_centroid(1,i-1,j) - abp_centroid(1,i,j) + x4 = abp_centroid(1,i+1,j) - abp_centroid(1,i,j) + x5 = abp_centroid(1,i+2,j) - abp_centroid(1,i,j) + + !IF (i == 1) THEN + ! x1 = piq + !ELSEIF (i == ncube_reconstruct-1) THEN + ! x5 = -piq + !ENDIF + + y1 = fcubehalo(i-2,j,p) + y2 = fcubehalo(i-1,j,p) + y3 = fcubehalo(i,j,p) + y4 = fcubehalo(i+1,j,p) + y5 = fcubehalo(i+2,j,p) + + denom(1) = (x2 - x1) * (x4 - x1) * (x5 - x1) * x1 + denom(2) = (x1 - x2) * (x4 - x2) * (x5 - x2) * x2 + denom(4) = (x1 - x4) * (x2 - x4) * (x5 - x4) * x4 + denom(5) = (x1 - x5) * (x2 - x5) * (x4 - x5) * x5 + + t(1) = x5 * x4 * x2 + t(2) = x5 * x4 * x1 + t(4) = x5 * x2 * x1 + t(5) = x4 * x2 * x1 + t(3) = (t(1) + t(2) + t(4) + t(5)) / (x1 * x2 * x4 * x5) + + pa(1) = x2 * x4 + x2 * x5 + x4 * x5 + pa(2) = x1 * x4 + x1 * x5 + x4 * x5 + pa(4) = x1 * x2 + x1 * x5 + x2 * x5 + pa(5) = x1 * x2 + x1 * x4 + x2 * x4 + pa(3) = (pa(1) + pa(2) + pa(4) + pa(5)) / (2.0 * x1 * x2 * x4 * x5) + + recons(1,i,j,p) = & + + y1 * t(1) / denom(1) & + + y2 * t(2) / denom(2) & + - y3 * t(3) & + + y4 * t(4) / denom(4) & + + y5 * t(5) / denom(5) + + IF (order > 2) THEN + recons(3,i,j,p) = & + - y1 * pa(1) / denom(1) & + - y2 * pa(2) / denom(2) & + + y3 * pa(3) & + - y4 * pa(4) / denom(4) & + - y5 * pa(5) / denom(5) + ENDIF + + ! Y-direction reconstruction + x1 = abp_centroid(2,i,j-2) - abp_centroid(2,i,j) + x2 = abp_centroid(2,i,j-1) - abp_centroid(2,i,j) + x4 = abp_centroid(2,i,j+1) - abp_centroid(2,i,j) + x5 = abp_centroid(2,i,j+2) - abp_centroid(2,i,j) + + !IF (j == 1) THEN + ! x1 = piq + !ELSEIF (j == ncube_reconstruct-1) THEN + ! x5 = -piq + !ENDIF + + y1 = fcubehalo(i,j-2,p) + y2 = fcubehalo(i,j-1,p) + y3 = fcubehalo(i,j,p) + y4 = fcubehalo(i,j+1,p) + y5 = fcubehalo(i,j+2,p) + + denom(1) = (x2 - x1) * (x4 - x1) * (x5 - x1) * x1 + denom(2) = (x1 - x2) * (x4 - x2) * (x5 - x2) * x2 + denom(4) = (x1 - x4) * (x2 - x4) * (x5 - x4) * x4 + denom(5) = (x1 - x5) * (x2 - x5) * (x4 - x5) * x5 + + t(1) = x5 * x4 * x2 + t(2) = x5 * x4 * x1 + t(4) = x5 * x2 * x1 + t(5) = x4 * x2 * x1 + t(3) = (t(1) + t(2) + t(4) + t(5)) / (x1 * x2 * x4 * x5) + + pa(1) = x2 * x4 + x2 * x5 + x4 * x5 + pa(2) = x1 * x4 + x1 * x5 + x4 * x5 + pa(4) = x1 * x2 + x1 * x5 + x2 * x5 + pa(5) = x1 * x2 + x1 * x4 + x2 * x4 + pa(3) = (pa(1) + pa(2) + pa(4) + pa(5)) / (2.0 * x1 * x2 * x4 * x5) + + recons(2,i,j,p) = & + + y1 * t(1) / denom(1) & + + y2 * t(2) / denom(2) & + - y3 * t(3) & + + y4 * t(4) / denom(4) & + + y5 * t(5) / denom(5) + + IF (order > 2) THEN + recons(4,i,j,p) = & + - y1 * pa(1) / denom(1) & + - y2 * pa(2) / denom(2) & + + y3 * pa(3) & + - y4 * pa(4) / denom(4) & + - y5 * pa(5) / denom(5) + recons(5,i,j,p) = 0.0 + ENDIF + + ENDDO + ENDDO + IF (order > 2) THEN + DO j = 1, ncube_reconstruct-1 + DO i = 1, ncube_reconstruct-1 + x1 = abp_centroid(1,i-1,j+1) - abp_centroid(1,i,j+1) + x2 = abp_centroid(1,i+1,j+1) - abp_centroid(1,i,j+1) + + y2 = (+ fcubehalo(i-1,j+1,p) * x2**2 & + - fcubehalo(i+1,j+1,p) * x1**2 & + - fcubehalo(i,j+1,p) * (x2**2 - x1**2)) / & + (x2 * x1 * (x2 - x1)) + + x1 = abp_centroid(1,i-1,j-1) - abp_centroid(1,i,j-1) + x2 = abp_centroid(1,i+1,j-1) - abp_centroid(1,i,j-1) + + y1 = (+ fcubehalo(i-1,j-1,p) * x2**2 & + - fcubehalo(i+1,j-1,p) * x1**2 & + - fcubehalo(i,j-1,p) * (x2**2 - x1**2)) / & + (x2 * x1 * (x2 - x1)) + + x1 = abp_centroid(2,i,j-1) - abp_centroid(2,i,j) + x2 = abp_centroid(2,i,j+1) - abp_centroid(2,i,j) + + recons(5,i,j,p) = & + (+ y1 * x2**2 & + - y2 * x1**2 & + - recons(1,i,j,p) * (x2**2 - x1**2)) / & + (x2 * x1 * (x2 - x1)) + + ENDDO + ENDDO + ENDIF + ENDDO + + END SUBROUTINE + +!------------------------------------------------------------------------------ +! SUBROUTINE ReconstructABPGradient_PLM +! +! Description: +! Construct a piecewise linear reconstruction of the gradient within +! each element on an ABP grid. +! +! Parameters: +! fcubehalo - Scalar field on the ABP grid to use in reconstruction +! recons (OUT) - Array of reconstructed coefficients for total elements +! order - Order of the scheme (2 or 3) +!------------------------------------------------------------------------------ + SUBROUTINE ReconstructABPGradient_PLM(fcubehalo, recons, order) + +! USE CubedSphereTrans +! USE InterpolateCSLL_Utils + + IMPLICIT NONE + + REAL (KIND=dbl_kind), & + DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6), INTENT(IN) :: fcubehalo + + REAL (KIND=dbl_kind), DIMENSION(:,:,:,:), INTENT(OUT) :: recons + + INTEGER (KIND=int_kind), INTENT(IN) :: order + + ! Local variables + INTEGER (KIND=int_kind) :: i, j, p + + REAL (KIND=dbl_kind) :: width + + ! ABP width between elements + width = pih / DBLE(ncube_reconstruct-1) + + DO p = 1, 6 + DO j = 1, ncube_reconstruct-1 + DO i = 1, ncube_reconstruct-1 + ! df/dx + recons(1,i,j,p) = (fcubehalo(i+1,j,p) - fcubehalo(i-1,j,p)) / & + (2.0 * width) + + ! df/dy + recons(2,i,j,p) = (fcubehalo(i,j+1,p) - fcubehalo(i,j-1,p)) / & + (2.0 * width) + + ! Stretching + recons(1,i,j,p) = recons(1,i,j,p) / (one + abp_centroid(1,i,j)**2) + recons(2,i,j,p) = recons(2,i,j,p) / (one + abp_centroid(2,i,j)**2) + + ! Third order scheme + IF (order > 2) THEN + ! d^2f/dx^2 + recons(3,i,j,p) = & + (fcubehalo(i+1,j,p) - 2.0 * fcubehalo(i,j,p) & + + fcubehalo(i-1,j,p)) / (width * width) + + ! d^2f/dy^2 + recons(4,i,j,p) = & + (fcubehalo(i,j+1,p) - 2.0 * fcubehalo(i,j,p) & + + fcubehalo(i,j-1,p)) / (width * width) + + ! d^2f/dxdy + recons(5,i,j,p) = & + (+ fcubehalo(i+1,j+1,p) - fcubehalo(i-1,j+1,p) & + - fcubehalo(i+1,j-1,p) + fcubehalo(i-1,j-1,p) & + ) / (4.0 * width * width) + + ! Stretching + recons(3,i,j,p) = & + (- 2.0 * abp_centroid(1,i,j) * (one + abp_centroid(1,i,j)**2) * recons(1,i,j,p) & + + recons(3,i,j,p)) / (one + abp_centroid(1,i,j)**2)**2 + + recons(4,i,j,p) = & + (- 2.0 * abp_centroid(2,i,j) * (one + abp_centroid(2,i,j)**2) * recons(2,i,j,p) & + + recons(4,i,j,p)) / (one + abp_centroid(2,i,j)**2)**2 + + recons(5,i,j,p) = recons(5,i,j,p) / & + ((one + abp_centroid(1,i,j)**2) * (one + abp_centroid(2,i,j)**2)) + + ! Scaling + recons(3,i,j,p) = 0.5 * recons(3,i,j,p) + recons(4,i,j,p) = 0.5 * recons(4,i,j,p) + + ENDIF + ENDDO + ENDDO + ENDDO + + END SUBROUTINE + +!------------------------------------------------------------------------------ +! SUBROUTINE ReconstructABPGradient_PPM +! +! Description: +! Construct a piecewise parabolic reconstruction of the gradient within +! each element on an ABP grid. +! +! Parameters: +! fcubehalo - Scalar field on the ABP grid to use in reconstruction +! recons (OUT) - Array of reconstructed coefficients for total elements +! order - Order of the scheme (2 or 3) +!------------------------------------------------------------------------------ + SUBROUTINE ReconstructABPGradient_PPM(fcubehalo, recons, order) + + +! USE CubedSphereTrans +! USE InterpolateCSLL_Utils + + IMPLICIT NONE + + REAL (KIND=dbl_kind), & + DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6), INTENT(IN) :: fcubehalo + + REAL (KIND=dbl_kind), DIMENSION(:,:,:,:), INTENT(OUT) :: recons + + INTEGER (KIND=int_kind), INTENT(IN) :: order + + ! Local variables + INTEGER (KIND=int_kind) :: i, j, p + + REAL (KIND=dbl_kind) :: width + + ! ABP width between elements + width = pih / DBLE(ncube_reconstruct-1) + + DO p = 1, 6 + DO j = 1, ncube_reconstruct-1 + DO i = 1, ncube_reconstruct-1 + ! df/dalfa + recons(1,i,j,p) = & + (+ fcubehalo(i+2,j,p) - 8.0 * fcubehalo(i+1,j,p) & + + 8.0 * fcubehalo(i-1,j,p) - fcubehalo(i-2,j,p)) / & + (- 12.0 * width) + + ! df/dbeta + recons(2,i,j,p) = & + (+ fcubehalo(i,j+2,p) - 8.0 * fcubehalo(i,j+1,p) & + + 8.0 * fcubehalo(i,j-1,p) - fcubehalo(i,j-2,p)) / & + (- 12.0 * width) + + ! Stretching + recons(1,i,j,p) = recons(1,i,j,p) / (one + abp_centroid(1,i,j)**2) + recons(2,i,j,p) = recons(2,i,j,p) / (one + abp_centroid(2,i,j)**2) + + ! Third order scheme + IF (order > 2) THEN + ! d^2f/dx^2 + recons(3,i,j,p) = (- fcubehalo(i+2,j,p) & + + 16_dbl_kind * fcubehalo(i+1,j,p) & + - 30_dbl_kind * fcubehalo(i,j,p) & + + 16_dbl_kind * fcubehalo(i-1,j,p) & + - fcubehalo(i-2,j,p) & + ) / (12_dbl_kind * width**2) + + ! d^2f/dy^2 + recons(4,i,j,p) = (- fcubehalo(i,j+2,p) & + + 16_dbl_kind * fcubehalo(i,j+1,p) & + - 30_dbl_kind * fcubehalo(i,j,p) & + + 16_dbl_kind * fcubehalo(i,j-1,p) & + - fcubehalo(i,j-2,p) & + ) / (12_dbl_kind * width**2) + + ! d^2f/dxdy + recons(5,i,j,p) = & + (+ fcubehalo(i+1,j+1,p) - fcubehalo(i-1,j+1,p) & + - fcubehalo(i+1,j-1,p) + fcubehalo(i-1,j-1,p) & + ) / (4.0 * width * width) + + ! Stretching + recons(3,i,j,p) = & + (- 2.0 * abp_centroid(1,i,j) * (one + abp_centroid(1,i,j)**2) * recons(1,i,j,p) & + + recons(3,i,j,p)) / (one + abp_centroid(1,i,j)**2)**2 + + recons(4,i,j,p) = & + (- 2.0 * abp_centroid(2,i,j) * (one + abp_centroid(2,i,j)**2) * recons(2,i,j,p) & + + recons(4,i,j,p)) / (one + abp_centroid(2,i,j)**2)**2 + + recons(5,i,j,p) = recons(5,i,j,p) / & + ((one + abp_centroid(1,i,j)**2) * (one + abp_centroid(2,i,j)**2)) + + ! Scaling + recons(3,i,j,p) = 0.5 * recons(3,i,j,p) + recons(4,i,j,p) = 0.5 * recons(4,i,j,p) + ENDIF + ENDDO + ENDDO + ENDDO + END SUBROUTINE + +!------------------------------------------------------------------------------ +! SUBROUTINE ReconstructABPGradient +! +! Description: +! Compute the reconstructed gradient in gnomonic coordinates for each +! ABP element. +! +! Parameters: +! fcube - Scalar field on the cubed sphere to use in reconstruction +! halomethod - Method for computing halo elements +! (0) Piecewise constant +! (1) Piecewise linear +! (3) Piecewise cubic +! recons_method - Method for computing the sub-grid scale gradient +! (0) Non-equidistant linear reconstruction +! (1) Non-equidistant parabolic reconstruction +! (2) Piecewise linear reconstruction with stretching +! (3) Piecewise parabolic reconstruction with stretching +! order - Order of the method being applied +! kmono - Apply monotone limiting (1) or not (0) +! recons (INOUT) - Array of reconstructed coefficients +!------------------------------------------------------------------------------ + SUBROUTINE ReconstructABPGradient( & + fcube, halomethod, recons_method, order, kmono, recons, kpd, kscheme) + +! USE InterpolateCSLL_Utils + + IMPLICIT NONE + + REAL (KIND=dbl_kind), & + DIMENSION(1:ncube_reconstruct-1, 1:ncube_reconstruct-1, 6), INTENT(IN) :: fcube + + INTEGER (KIND=int_kind), INTENT(IN) :: halomethod, recons_method + INTEGER (KIND=int_kind), INTENT(IN) :: order, kmono, kpd, kscheme + + REAL (KIND=dbl_kind), DIMENSION(:,:,:,:), INTENT(INOUT) :: recons + + ! Local variables + INTEGER (KIND=int_kind) :: i, j, p + + REAL (KIND=dbl_kind), DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6) :: fcubehalo + + ! Report status + WRITE (*,*) '...Performing sub-grid scale reconstruction on ABP grid' + + ! Compute element haloes + WRITE(*,*) "fill cubed-sphere halo for reconstruction" + DO p = 1, 6 + IF (halomethod == 0) THEN + CALL CubedSphereFillHalo(fcube, fcubehalo, p, ncube_reconstruct, 2) + + ELSEIF (halomethod == 1) THEN + CALL CubedSphereFillHalo_Linear(fcube, fcubehalo, p, ncube_reconstruct) + + ELSEIF (halomethod == 3) THEN + !halomethod is always 3 in the standard CSLAM setup + CALL CubedSphereFillHalo_Cubic(fcube, fcubehalo, p, ncube_reconstruct) + ELSE + WRITE (*,*) 'Fatal Error: In ReconstructABPGradient' + WRITE (*,*) 'Invalid halo method: ', halomethod + WRITE (*,*) 'Halo method must be 0, 1 or 3.' + STOP + ENDIF + ENDDO + + ! Nonequidistant linear reconstruction + IF (recons_method == 1) THEN + CALL ReconstructABPGradient_NEL(fcubehalo, recons, order) + + ! Nonequidistant parabolic reconstruction (JCP paper) + ELSEIF (recons_method == 2) THEN + WRITE(*,*) "Nonequidistant parabolic reconstruction" + CALL ReconstructABPGradient_NEP(fcubehalo, recons, order) + + ! Piecewise linear reconstruction with rotation + ELSEIF (recons_method == 3) THEN + CALL ReconstructABPGradient_PLM(fcubehalo, recons, order) + + ! Piecewise parabolic reconstruction with rotation + ELSEIF (recons_method == 4) THEN + CALL ReconstructABPGradient_PPM(fcubehalo, recons, order) + + ELSE + WRITE(*,*) 'Fatal Error: In ReconstructABPGradient' + WRITE(*,*) 'Specified recons_method out of range. Given: ', recons_method + WRITE(*,*) 'Valid values: 1, 2, 3, 4' + STOP + ENDIF + + ! Apply monotone filtering + SELECT CASE (kmono) + CASE (0) !Do nothing + WRITE(*,*) "no filter applied to the reconstruction" + CASE (1) + + !Simplest filter: just scales the recon so it's extreme value + !is no bigger than the original values of this point and its neighbors + CALL MonotonizeABPGradient(fcubehalo, order, recons, .FALSE.) + + CASE (2) + + !Applies a more sophisticated Van Leer limiter (or, to be consistent, a filter) + CALL VanLeerLimit(fcubehalo, order, recons) + + CASE (3) + + !Applies a selective filter + CALL MonotonizeABPGradient(fcubehalo, order, recons, .TRUE.) + + CASE (4) + + !A filter that filters the linear part first + CALL MonotonizeABPGradient_New(fcubehalo, order, recons) + + CASE DEFAULT + WRITE(*,*) "Limiter kmono = ", kmono, " does not exist." + STOP 1201 + + END SELECT + + !Apply positive-definite filtering, if desired. This should + !ONLY be applied to the S-L method, since the flux-form + !method needs something different done. (In particular, using + !positive-definite reconstructions does not ensure that a flux- + !form scheme is positive definite, since we could get negatives + !when subtracting the resulting fluxes.) + !HOWEVER...we will allow this to be enabled, for testing purposes + IF ( (kpd > 0 .AND. kscheme == 2) .OR. (kpd == 2 .AND. kscheme == 4) ) THEN + WRITE(*,*) "applying positive deifnite constraint" + CALL PosDefABPGradient(fcubehalo, order, recons) + END IF + + + END SUBROUTINE + + + +!------------------------------------------------------------------------------ +!------------------------------------------------------------------------------ +! SUBROUTINE AdjustLimiter +! +! Description: +! Adjust the slope limiter based on new point values. +! +! Parameters: +! value - Point value +! element_value - Value at the center of the element +! local_max - Local maximum value of the function (from neighbours) +! local_min - Local minimum value of the function (to neighbours) +! min_phi (INOUT) - Slope limiter +!------------------------------------------------------------------------------ + SUBROUTINE AdjustLimiter(value, element_value, local_min, local_max, min_phi) + + IMPLICIT NONE + + REAL (KIND=dbl_kind), INTENT(IN) :: value, element_value + REAL (KIND=dbl_kind), INTENT(IN) :: local_min, local_max + REAL (KIND=dbl_kind), INTENT(INOUT) :: min_phi + + ! Local variables + REAL (KIND=dbl_kind) :: phi = 0.0 + + IF ((local_min > element_value ) .OR. (local_max < element_value )) THEN + WRITE (*,*) 'Fatal Error: In AdjustLimiter' + WRITE (*,*) 'Local min: ', local_min, ' max: ', local_max + WRITE (*,*) 'Elemn: ', element_value + STOP + ENDIF + + ! Check against the minimum bound on the reconstruction + IF (value - element_value > tiny * value) THEN + phi = (local_max - element_value) / & + (value - element_value) + + min_phi = MIN(min_phi, phi) + + ! Check against the maximum bound on the reconstruction + ELSEIF (value - element_value < -tiny * value) THEN + phi = (local_min - element_value) / & + (value - element_value) + + min_phi = MIN(min_phi, phi) + + ENDIF + + IF (min_phi < 0.0) THEN + WRITE (*,*) 'Fatal Error: In AdjustLimiter' + WRITE (*,*) 'Min_Phi: ', min_phi + WRITE (*,*) 'Phi: ', phi + WRITE (*,*) 'Value: ', value + WRITE (*,*) 'Elemn: ', element_value + WRITE (*,*) 'Val-E: ', value - element_value + STOP + ENDIF + + END SUBROUTINE + +!------------------------------------------------------------------------------ +! SUBROUTINE VanLeerLimit +! +! Description: +! Apply a 2D Van Leer-type limiter to a reconstruction. This acts ONLY +! on the linear part of the reconstruction , if any. If passed a PCoM +! reconstruction, this just returns without altering the recon. +! +! Parameters: +! fcubehalo - Scalar field on the cubed sphere to use in reconstruction +! order - Order of the reconstruction +! recons (INOUT) - Array of reconstructed coefficients +! +! Remarks: +! The Van Leer Limiter described here is given on pages 328--329 +! of Dukowicz and Baumgardner (2000). There are no guarantees +! on what it will do to PPM. +!------------------------------------------------------------------------------ + SUBROUTINE VanLeerLimit(fcubehalo, order, recons) + + + IMPLICIT NONE + + REAL (KIND=dbl_kind), DIMENSION(-1:ncube_reconstruct+1, -1:ncube_reconstruct+1, 6), & + INTENT(IN) :: fcubehalo + + INTEGER (KIND=int_kind), INTENT(IN) :: order + + REAL (KIND=dbl_kind), DIMENSION(:,:,:,:), INTENT(INOUT) :: recons + + ! Local variables + INTEGER (KIND=int_kind) :: i, j, k, m, n + + REAL (KIND=dbl_kind) :: local_min, local_max, value, phi, min_phi, & + recon_min, recon_max + + ! The first-order piecewise constant scheme is monotone by construction + IF (order == 1) THEN + RETURN + ENDIF + + ! Apply monotone limiting + DO k = 1, 6 + DO j = 1, ncube_reconstruct-1 + DO i = 1, ncube_reconstruct-1 + CALL ABPHaloMinMax(fcubehalo, i, j, k, local_min, local_max,.FALSE.) + + ! Initialize the limiter + min_phi = one + + ! For the second-order calculation, the minima and maxima will occur + ! at the corner points of the element. For the Van Leer limiter, we + !wish to find BOTH of the reconstruction extrema. + recon_min = bignum + recon_max = -bignum + + DO m = i, i+1 + DO n = j, j+1 + + ! Evaluate the function at each corner point + CALL EvaluateABPReconstruction( & + fcubehalo, recons, i, j, k, gp(m), gp(n), order, value) + recon_min = MIN(recon_min, value) + recon_max = MAX(recon_max, value) + + ENDDO + ENDDO + + !This is equation 27 in Dukowicz and Baumgardner 2000 + min_phi = MIN(one, MAX(0.0, (local_min - fcubehalo(i,j,k))/(recon_min - fcubehalo(i,j,k))), & + MAX(0.0, (local_max - fcubehalo(i,j,k))/(recon_max - fcubehalo(i,j,k))) ) + + IF ((min_phi < -tiny) .OR. (min_phi > one + tiny)) THEN + WRITE (*,*) 'Fatal Error: In MonotonizeABPGradient' + WRITE (*,*) 'Slope limiter out of range: ', min_phi + STOP + ENDIF + + ! Apply monotone limiter to all reconstruction coefficients + recons(1,i,j,k) = min_phi * recons(1,i,j,k) + recons(2,i,j,k) = min_phi * recons(2,i,j,k) + + END DO + END DO + END DO + + + + + END SUBROUTINE VanLeerLimit + + !------------------------------------------------------------------------------ + ! SUBROUTINE EquiangularElementArea + ! + ! Description: + ! Compute the area of a single equiangular cubed sphere grid cell. + ! + ! Parameters: + ! alpha - Alpha coordinate of lower-left corner of grid cell + ! da - Delta alpha + ! beta - Beta coordinate of lower-left corner of grid cell + ! db - Delta beta + !------------------------------------------------------------------------------ + REAL(KIND=dbl_kind) FUNCTION EquiangularElementArea(alpha, da, beta, db) + + IMPLICIT NONE + +! REAL (kind=dbl_kind) :: EquiangularElementArea + REAL (kind=dbl_kind) :: alpha, da, beta, db + REAL (kind=dbl_kind) :: a1, a2, a3, a4 + + ! Calculate interior grid angles + a1 = EquiangularGridAngle(alpha , beta ) + a2 = pi - EquiangularGridAngle(alpha+da, beta ) + a3 = pi - EquiangularGridAngle(alpha , beta+db) + a4 = EquiangularGridAngle(alpha+da, beta+db) + + ! Area = r*r*(-2*pi+sum(interior angles)) + EquiangularElementArea = -pi2 + a1 + a2 + a3 + a4 + + END FUNCTION EquiangularElementArea + + !------------------------------------------------------------------------------ + ! FUNCTION EquiangularGridAngle + ! + ! Description: + ! Compute the angle between equiangular cubed sphere projection grid lines. + ! + ! Parameters: + ! alpha - Alpha coordinate of evaluation point + ! beta - Beta coordinate of evaluation point + !------------------------------------------------------------------------------ + REAL(KIND=dbl_kind) FUNCTION EquiangularGridAngle(alpha, beta) + IMPLICIT NONE + REAL (kind=dbl_kind) :: alpha, beta + EquiangularGridAngle = ACOS(-SIN(alpha) * SIN(beta)) + END FUNCTION EquiangularGridAngle + +!------------------------------------------------------------------------------ +! SUBROUTINE CubedSphereFillHalo +! +! Description: +! Recompute the cubed sphere data storage array, with the addition of a +! halo region around the specified panel. +! +! Parameters: +! parg - Current panel values +! zarg (OUT) - Calculated panel values with halo/ghost region +! np - Panel number +! ncube - Dimension of the cubed sphere (# of grid lines) +! nhalo - Number of halo/ghost elements around each panel +!------------------------------------------------------------------------------ + SUBROUTINE CubedSphereFillHalo(parg, zarg, np, ncube, nhalo) + + IMPLICIT NONE + + REAL (KIND=dbl_kind), DIMENSION(ncube-1, ncube-1, 6), INTENT(IN) :: parg + + REAL (KIND=dbl_kind), & + DIMENSION(1-nhalo:ncube+nhalo-1, 1-nhalo:ncube+nhalo-1, 6), & + INTENT(OUT) :: zarg + + INTEGER (KIND=int_kind), INTENT(IN) :: np, ncube,nhalo + + ! Local variables + INTEGER (KIND=int_kind) :: jh,jhy + + !zarg = 0.0 !DBG + zarg(1:ncube-1,1:ncube-1,np) = parg(1:ncube-1,1:ncube-1,np) + + zarg(1-nhalo:0,1-nhalo:0,np) = 0.0 + zarg(1-nhalo:0,ncube:ncube+nhalo-1,np) = 0.0 + zarg(ncube:ncube+nhalo-1,1-nhalo:0,np) = 0.0 + zarg(ncube:ncube+nhalo-1,ncube:ncube+nhalo-1,np) = 0.0 + + ! Equatorial panels + IF (np==1) THEN + DO jh=1,nhalo + zarg(ncube+jh-1,1:ncube-1 ,1) = parg(jh ,1:ncube-1 ,2) !exchange right + zarg(1-jh ,1:ncube-1 ,1) = parg(ncube-jh ,1:ncube-1 ,4) !exchange left + zarg(1:ncube-1 ,1-jh ,1) = parg(1:ncube-1 ,ncube-jh ,5) !exchange below + zarg(1:ncube-1 ,ncube+jh-1,1) = parg(1:ncube-1 ,jh ,6) !exchange over + ENDDO + + ELSE IF (np==2) THEN + DO jh=1,nhalo + zarg(1-jh ,1:ncube-1 ,2) = parg(ncube-jh,1:ncube-1 ,1) !exchange left + zarg(ncube+jh-1,1:ncube-1 ,2) = parg(jh ,1:ncube-1 ,3) !exchange right + zarg(1:ncube-1 ,1-jh ,2) = parg(ncube-jh,ncube-1:1:-1,5) !exchange below + zarg(1:ncube-1 ,ncube+jh-1,2) = parg(ncube-jh,1:ncube-1 ,6) !exchange over + ENDDO + + ELSE IF (np==3) THEN + DO jh=1,nhalo + zarg(ncube+jh-1,1:ncube-1 ,3) = parg(jh ,1:ncube-1,4) !exchange right + zarg(1-jh ,1:ncube-1 ,3) = parg(ncube-jh ,1:ncube-1,2) !exchange left + zarg(1:ncube-1 ,1-jh ,3) = parg(ncube-1:1:-1,jh ,5) !exchange below + zarg(1:ncube-1 ,ncube+jh-1,3) = parg(ncube-1:1:-1,ncube-jh ,6) !exchange over + ENDDO + + ELSE IF (np==4) THEN + DO jh=1,nhalo + zarg(1-jh ,1:ncube-1 ,4) = parg(ncube-jh,1:ncube-1 ,3) !exchange left + zarg(ncube+jh-1,1:ncube-1 ,4) = parg(jh ,1:ncube-1 ,1) !exchange right + zarg(1:ncube-1 ,1-jh ,4) = parg(jh ,1:ncube-1 ,5) !exchange below + zarg(1:ncube-1 ,ncube+jh-1,4) = parg(jh ,ncube-1:1:-1,6) !exchange over + ENDDO + + ! Bottom panel + ELSE IF (np==5) THEN + DO jh=1,nhalo + zarg(1-jh ,1:ncube-1 ,5) = parg(1:ncube-1 ,jh ,4) !exchange left + zarg(ncube+jh-1,1:ncube-1 ,5) = parg(ncube-1:1:-1,jh ,2) !exchange right + zarg(1:ncube-1 ,1-jh ,5) = parg(ncube-1:1:-1,jh ,3) !exchange below + zarg(1:ncube-1 ,ncube+jh-1,5) = parg(1:ncube-1 ,jh ,1) !exchange over + ENDDO + + ! Top panel + ELSE IF (np==6) THEN + DO jh=1,nhalo + zarg(1-jh ,1:ncube-1 ,6) = parg(ncube-1:1:-1,ncube-jh,4) !exchange left + zarg(ncube+jh-1,1:ncube-1 ,6) = parg(1:ncube-1 ,ncube-jh,2) !exchange right + zarg(1:ncube-1 ,1-jh ,6) = parg(1:ncube-1 ,ncube-jh,1) !exchange below + zarg(1:ncube-1 ,ncube+jh-1,6) = parg(ncube-1:1:-1,ncube-jh,3) !exchange over + ENDDO + + ELSE + WRITE (*,*) 'Fatal error: In CubedSphereFillHalo' + WRITE (*,*) 'Invalid panel id ', np + STOP + ENDIF + + END SUBROUTINE CubedSphereFillHalo + +!------------------------------------------------------------------------------ +! SUBROUTINE CubedSphereFillHalo_Linear +! +! Description: +! Recompute the cubed sphere data storage array, with the addition of a +! 2-element halo region around the specified panel. Use linear order +! interpolation to translate between panels. +! +! Parameters: +! parg - Current panel values +! zarg (OUT) - Calculated panel values with halo/ghost region +! np - Panel number +! ncube - Dimension of the cubed sphere (# of grid lines) +!------------------------------------------------------------------------------ + SUBROUTINE CubedSphereFillHalo_Linear(parg, zarg, np, ncube) + +! USE CubedSphereTrans ! Cubed sphere transforms + + IMPLICIT NONE + + INTEGER (KIND=int_kind), PARAMETER :: nhalo = 2 + + REAL (KIND=dbl_kind), DIMENSION(ncube-1, ncube-1, 6), INTENT(IN) :: parg + + REAL (KIND=dbl_kind), & + DIMENSION(1-nhalo:ncube+nhalo-1, 1-nhalo:ncube+nhalo-1, 6), & + INTENT(OUT) :: zarg + + INTEGER (KIND=int_kind), INTENT(IN) :: np, ncube + + ! Local variables + INTEGER (KIND=int_kind) :: ii, iref, jj, ipanel, imin, imax + REAL (KIND=dbl_kind) :: width, lon, lat, beta, a, newbeta + + REAL (KIND=dbl_kind), DIMENSION(0:ncube, nhalo) :: prealpha + REAL (KIND=dbl_kind), DIMENSION(0:ncube, nhalo) :: newalpha + + REAL (KIND=dbl_kind), & + DIMENSION(1-nhalo:ncube+nhalo-1, 1-nhalo:ncube+nhalo-1, 6) :: yarg + + ! Use 0.0 order interpolation to begin + CALL CubedSphereFillHalo(parg, yarg, np, ncube, nhalo) + + zarg(:,:,np) = yarg(:,:,np) + + ! Calculate the overlapping alpha coordinates + width = pih / DBLE(ncube-1) + + DO jj = 1, nhalo + DO ii = 0, ncube + prealpha(ii, jj) = width * (DBLE(ii-1) + 0.5) - piq + beta = - width * (DBLE(jj-1) + 0.5) - piq + + CALL CubedSphereABPFromABP(prealpha(ii,jj), beta, 1, 5, & + newalpha(ii,jj), newbeta) + ENDDO + ENDDO + + ! Now apply linear interpolation to obtain edge components + DO jj = 1, nhalo + ! Reset the reference index + iref = 2 + + ! Interpolation can be applied to more elements after first band + IF (jj == 1) THEN + imin = 1 + imax = ncube-1 + ELSE + imin = 0 + imax = ncube + ENDIF + + ! Apply linear interpolation + DO ii = imin, imax + DO WHILE ((iref .NE. ncube-1) .AND. & + (newalpha(ii,jj) > prealpha(iref,jj))) + iref = iref + 1 + ENDDO + + IF ((newalpha(ii,jj) > prealpha(iref-1,jj)) .AND. & + (newalpha(ii,jj) .LE. prealpha(iref ,jj))) & + THEN + a = (newalpha(ii,jj) - prealpha(iref-1,jj)) / & + (prealpha(iref,jj) - prealpha(iref-1,jj)) + + IF ((a < 0.0) .OR. (a > one)) THEN + WRITE (*,*) 'FAIL in CubedSphereFillHalo_Linear' + WRITE (*,*) 'a out of bounds' + STOP + ENDIF + + ! Bottom edge of panel + zarg(ii, 1-jj, np) = & + (one - a) * yarg(iref-1, 1-jj, np) + & + a * yarg(iref, 1-jj, np) + + ! Left edge of panel + zarg(1-jj, ii, np) = & + (one - a) * yarg(1-jj, iref-1, np) + & + a * yarg(1-jj, iref, np) + + ! Top edge of panel + zarg(ii, ncube+jj-1, np) = & + (one - a) * yarg(iref-1, ncube+jj-1, np) + & + a * yarg(iref, ncube+jj-1, np) + + ! Right edge of panel + zarg(ncube+jj-1, ii, np) = & + (one - a) * yarg(ncube+jj-1, iref-1, np) + & + a * yarg(ncube+jj-1, iref, np) + + ELSE + WRITE (*,*) 'FAIL in CubedSphereFillHalo_Linear' + WRITE (*,*) 'ii: ', ii, ' jj: ', jj + WRITE (*,*) 'newalpha: ', newalpha(ii,jj) + WRITE (*,*) 'prealpha: ', prealpha(iref-1,jj), '-', prealpha(iref,jj) + STOP + ENDIF + ENDDO + ENDDO + + ! Fill in corner bits + zarg(0, 0, np) = & + 0.25 * (zarg(1,0,np) + zarg(0,1,np) + & + zarg(-1,0,np) + zarg(0,-1,np)) + zarg(0, ncube, np) = & + 0.25 * (zarg(0,ncube-1,np) + zarg(0,ncube+1,np) + & + zarg(-1,ncube,np) + zarg(1,ncube,np)) + zarg(ncube, 0, np) = & + 0.25 * (zarg(ncube-1,0,np) + zarg(ncube+1,0,np) + & + zarg(ncube,-1,np) + zarg(ncube,1,np)) + zarg(ncube, ncube, np) = & + 0.25 * (zarg(ncube-1,ncube,np) + zarg(ncube+1,ncube,np) + & + zarg(ncube,ncube-1,np) + zarg(ncube,ncube+1,np)) + + END SUBROUTINE CubedSphereFillHalo_Linear + +!------------------------------------------------------------------------------ +! SUBROUTINE CubedSphereFillHalo_Cubic +! +! Description: +! Recompute the cubed sphere data storage array, with the addition of a +! 2-element halo region around the specified panel. Use higher order +! interpolation to translate between panels. +! +! Parameters: +! parg - Current panel values +! zarg (OUT) - Calculated panel values with halo/ghost region +! np - Panel number +! ncube - Dimension of the cubed sphere (# of grid lines) +!------------------------------------------------------------------------------ + SUBROUTINE CubedSphereFillHalo_Cubic(parg, zarg, np, ncube) + +! USE CubedSphereTrans ! Cubed sphere transforms +! USE MathUtils ! Has function for 1D cubic interpolation + + IMPLICIT NONE + + INTEGER (KIND=int_kind), PARAMETER :: nhalo = 2 + + REAL (KIND=dbl_kind), DIMENSION(ncube-1, ncube-1, 6), INTENT(IN) :: parg + + REAL (KIND=dbl_kind), & + DIMENSION(1-nhalo:ncube+nhalo-1, 1-nhalo:ncube+nhalo-1, 6), & + INTENT(OUT) :: zarg + + INTEGER (KIND=int_kind), INTENT(IN) :: np, ncube + + ! Local variables + INTEGER (KIND=int_kind) :: ii, iref, ibaseref, jj, ipanel, imin, imax + REAL (KIND=dbl_kind) :: width, lon, lat, beta, a, newbeta + + REAL (KIND=dbl_kind), DIMENSION(0:ncube, nhalo) :: prealpha + REAL (KIND=dbl_kind), DIMENSION(0:ncube, nhalo) :: newalpha + REAL (KIND=dbl_kind), DIMENSION(1:4) :: C, D, X + + REAL (KIND=dbl_kind), & + DIMENSION(1-nhalo:ncube+nhalo-1, 1-nhalo:ncube+nhalo-1, 6) :: yarg + + ! Use 0.0 order interpolation to begin + CALL CubedSphereFillHalo(parg, yarg, np, ncube, nhalo) + + zarg(:,:,np) = yarg(:,:,np) + + ! Calculate the overlapping alpha coordinates + width = pih / DBLE(ncube-1) + + DO jj = 1, nhalo + DO ii = 0, ncube + ! + ! alpha,beta for the cell center (extending the panel) + ! + prealpha(ii, jj) = width * (DBLE(ii-1) + 0.5) - piq + beta = - width * (DBLE(jj-1) + 0.5) - piq + + CALL CubedSphereABPFromABP(prealpha(ii,jj), beta, 1, 5, & + newalpha(ii,jj), newbeta) + ENDDO + ENDDO + + ! Now apply cubic interpolation to obtain edge components + DO jj = 1, nhalo + ! Reset the reference index, which gives the element in newalpha that + ! is closest to ii, looking towards larger values of alpha. + iref = 2 + + ! Interpolation can be applied to more elements after first band +! IF (jj == 1) THEN +! imin = 1 +! imax = ncube-1 +! ELSE + imin = 0 + imax = ncube +! ENDIF + + ! Apply cubic interpolation + DO ii = imin, imax + DO WHILE ((iref .NE. ncube-1) .AND. & + (newalpha(ii,jj) > prealpha(iref,jj))) + iref = iref + 1 + ENDDO + + ! Smallest index for cubic interpolation - apply special consideration + IF (iref == 2) THEN + ibaseref = iref-1 + + ! Largest index for cubic interpolation - apply special consideration + ELSEIF (iref == ncube-1) THEN + ibaseref = iref-3 + + ! Normal range + ELSE + ibaseref = iref-2 + ENDIF + + ! Bottom edge of panel + zarg(ii, 1-jj, np) = & + CUBIC_EQUISPACE_INTERP( & + width, newalpha(ii,jj) - prealpha(ibaseref,jj), & + yarg(ibaseref:ibaseref+3, 1-jj, np)) + + ! Left edge of panel + zarg(1-jj, ii, np) = & + CUBIC_EQUISPACE_INTERP( & + width, newalpha(ii,jj) - prealpha(ibaseref,jj), & + yarg(1-jj, ibaseref:ibaseref+3, np)) + + ! Top edge of panel + zarg(ii, ncube+jj-1, np) = & + CUBIC_EQUISPACE_INTERP( & + width, newalpha(ii,jj) - prealpha(ibaseref,jj), & + yarg(ibaseref:ibaseref+3, ncube+jj-1, np)) + + ! Right edge of panel + zarg(ncube+jj-1, ii, np) = & + CUBIC_EQUISPACE_INTERP( & + width, newalpha(ii,jj) - prealpha(ibaseref,jj), & + yarg(ncube+jj-1, ibaseref:ibaseref+3, np)) + + ENDDO + ENDDO + + ! Fill in corner bits + zarg(0, 0, np) = & + 0.25 * (zarg(1,0,np) + zarg(0,1,np) + & + zarg(-1,0,np) + zarg(0,-1,np)) + zarg(0, ncube, np) = & + 0.25 * (zarg(0,ncube-1,np) + zarg(0,ncube+1,np) + & + zarg(-1,ncube,np) + zarg(1,ncube,np)) + zarg(ncube, 0, np) = & + 0.25 * (zarg(ncube-1,0,np) + zarg(ncube+1,0,np) + & + zarg(ncube,-1,np) + zarg(ncube,1,np)) + zarg(ncube, ncube, np) = & + 0.25 * (zarg(ncube-1,ncube,np) + zarg(ncube+1,ncube,np) + & + zarg(ncube,ncube-1,np) + zarg(ncube,ncube+1,np)) + + END SUBROUTINE CubedSphereFillHalo_Cubic + +!------------------------------------------------------------------------------ +! SUBROUTINE CubedSphereABPFromABP +! +! Description: +! Determine the (alpha,beta,idest) coordinate of a source point on +! panel isource. +! +! Parameters: +! alpha_in - Alpha coordinate in +! beta_in - Beta coordinate in +! isource - Source panel +! idest - Destination panel +! alpha_out (OUT) - Alpha coordinate out +! beta_out (OUT) - Beta coordiante out +!------------------------------------------------------------------------------ + SUBROUTINE CubedSphereABPFromABP(alpha_in, beta_in, isource, idest, & + alpha_out, beta_out) + + IMPLICIT NONE + + REAL (KIND=dbl_kind), INTENT(IN) :: alpha_in, beta_in + INTEGER (KIND=int_kind), INTENT(IN) :: isource, idest + REAL (KIND=dbl_kind), INTENT(OUT) :: alpha_out, beta_out + + ! Local variables + REAL (KIND=dbl_kind) :: a1, b1 + REAL (KIND=dbl_kind) :: xx, yy, zz + REAL (KIND=dbl_kind) :: sx, sy, sz + + ! Convert to relative Cartesian coordinates + a1 = TAN(alpha_in) + b1 = TAN(beta_in) + + sz = (one + a1 * a1 + b1 * b1)**(-0.5) + sx = sz * a1 + sy = sz * b1 + + ! Convert to full Cartesian coordinates + IF (isource == 6) THEN + yy = sx; xx = -sy; zz = sz + + ELSEIF (isource == 5) THEN + yy = sx; xx = sy; zz = -sz + + ELSEIF (isource == 1) THEN + yy = sx; zz = sy; xx = sz + + ELSEIF (isource == 3) THEN + yy = -sx; zz = sy; xx = -sz + + ELSEIF (isource == 2) THEN + xx = -sx; zz = sy; yy = sz + + ELSEIF (isource == 4) THEN + xx = sx; zz = sy; yy = -sz + + ELSE + WRITE(*,*) 'Fatal Error: Source panel invalid in CubedSphereABPFromABP' + WRITE(*,*) 'panel = ', isource + STOP + ENDIF + + ! Convert to relative Cartesian coordinates on destination panel + IF (idest == 6) THEN + sx = yy; sy = -xx; sz = zz + + ELSEIF (idest == 5) THEN + sx = yy; sy = xx; sz = -zz + + ELSEIF (idest == 1) THEN + sx = yy; sy = zz; sz = xx + + ELSEIF (idest == 3) THEN + sx = -yy; sy = zz; sz = -xx + + ELSEIF (idest == 2) THEN + sx = -xx; sy = zz; sz = yy + + ELSEIF (idest == 4) THEN + sx = xx; sy = zz; sz = -yy + + ELSE + WRITE(*,*) 'Fatal Error: Dest panel invalid in CubedSphereABPFromABP' + WRITE(*,*) 'panel = ', idest + STOP + ENDIF + IF (sz < 0) THEN + WRITE(*,*) 'Fatal Error: In CubedSphereABPFromABP' + WRITE(*,*) 'Invalid relative Z coordinate' + STOP + ENDIF + + ! Use panel information to calculate (alpha, beta) coords + alpha_out = ATAN(sx / sz) + beta_out = ATAN(sy / sz) + + END SUBROUTINE + + +!------------------------------------------------------------------------------ +! FUNCTION CUBIC_EQUISPACE_INTERP +! +! Description: +! Apply cubic interpolation on the specified array of values, where all +! points are equally spaced. +! +! Parameters: +! dx - Spacing of points +! x - X coordinate where interpolation is to be applied +! y - Array of 4 values = f(x + k * dx) where k = 0,1,2,3 +!------------------------------------------------------------------------------ + FUNCTION CUBIC_EQUISPACE_INTERP(dx, x, y) + + IMPLICIT NONE + + REAL (KIND=dbl_kind) :: CUBIC_EQUISPACE_INTERP + REAL (KIND=dbl_kind) :: dx, x + REAL (KIND=dbl_kind), DIMENSION(1:4) :: y + + CUBIC_EQUISPACE_INTERP = & + (-y(1) / (6.0 * dx**3)) * (x - dx) * (x - 2.0 * dx) * (x - 3.0 * dx) + & + ( y(2) / (2.0 * dx**3)) * (x) * (x - 2.0 * dx) * (x - 3.0 * dx) + & + (-y(3) / (2.0 * dx**3)) * (x) * (x - dx) * (x - 3.0 * dx) + & + ( y(4) / (6.0 * dx**3)) * (x) * (x - dx) * (x - 2.0 * dx) + + END FUNCTION CUBIC_EQUISPACE_INTERP + +! FUNCTION I_10_ab(alpha,beta) +! IMPLICIT NONE +! REAL (KIND=dbl_kind) :: I_10_AB +! REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta +! I_10_ab = -ASINH(COS(alpha) * TAN(beta)) +! END FUNCTION I_10_AB +!! +! +! REAL (KIND=dbl_kind) FUNCTION I_01_ab(alpha,beta) +! IMPLICIT NONE +! REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta +! I_01_ab = -ASINH(COS(beta) * TAN(alpha)) +! END FUNCTION I_01_AB +! +! REAL (KIND=dbl_kind) FUNCTION I_20_ab(alpha,beta) +! IMPLICIT NONE +! REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta +! +! I_20_ab = TAN(beta)*ASINH(COS(beta)*TAN(alpha))+ACOS(SIN(alpha)*SIN(beta)) +! END FUNCTION I_20_AB +! +! REAL (KIND=dbl_kind) FUNCTION I_02_ab(alpha,beta) +! IMPLICIT NONE +! REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta +! +! I_02_ab = TAN(alpha)*ASINH(TAN(beta)*COS(alpha))+ACOS(SIN(alpha)*SIN(beta)) +! END FUNCTION I_02_AB +! +! REAL (KIND=dbl_kind) FUNCTION I_11_ab(alpha,beta) +! IMPLICIT NONE +! REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta +! +! I_11_ab = -SQRT(1.0+TAN(alpha)**2+TAN(beta)**2) +! END FUNCTION I_11_AB +! + + +END MODULE reconstruct + diff --git a/components/eam/tools/orographic_drag_toolkit/remap.F90 b/components/eam/tools/orographic_drag_toolkit/remap.F90 new file mode 100755 index 000000000000..fbfca57006a9 --- /dev/null +++ b/components/eam/tools/orographic_drag_toolkit/remap.F90 @@ -0,0 +1,1564 @@ +MODULE remap + INTEGER, PARAMETER :: & + int_kind = KIND(1), & + real_kind = SELECTED_REAL_KIND(p=14,r=100),& + dbl_kind = selected_real_kind(13) + + INTEGER :: nc,nhe + +! LOGICAL, PARAMETER:: ldbgr_r = .FALSE. + LOGICAL :: ldbgr + LOGICAL :: ldbg_global + + REAL(kind=real_kind), PARAMETER :: & + one = 1.0 ,& + aa = 1.0 ,& + tiny= 1.0E-9 ,& + bignum = 1.0E20 + REAL (KIND=dbl_kind), parameter :: fuzzy_width = 10.0*tiny !CAM-SE add + + contains + + + subroutine compute_weights_cell(xcell_in,ycell_in,jx,jy,nreconstruction,xgno,ygno,& + jx_min, jx_max, jy_min, jy_max,tmp,& + ngauss,gauss_weights,abscissae,weights,weights_eul_index,jcollect,jmax_segments,& + nc_in,nhe_in,nvertex,ldbg) + + implicit none + integer (kind=int_kind) , intent(in):: nreconstruction, jx,jy,ngauss,jmax_segments + real (kind=real_kind) , dimension(0:nvertex+1) :: xcell_in,ycell_in +! real (kind=real_kind) , dimension(0:5), intent(in):: xcell_in,ycell_in + integer (kind=int_kind), intent(in) :: nc_in,nhe_in,nvertex + logical, intent(in) :: ldbg + ! + ! ipanel is just for debugging + ! + integer (kind=int_kind), intent(in) :: jx_min, jy_min, jx_max, jy_max + real (kind=real_kind), dimension(-nhe_in:nc_in+2+nhe_in), intent(in) :: xgno + real (kind=real_kind), dimension(-nhe_in:nc_in+2+nhe_in), intent(in) :: ygno + ! + ! for Gaussian quadrature + ! + real (kind=real_kind), dimension(ngauss), intent(in) :: gauss_weights, abscissae + ! + ! boundaries of domain + ! + real (kind=real_kind):: tmp + ! + ! Number of Eulerian sub-cell integrals for the cell in question + ! + integer (kind=int_kind), intent(out) :: jcollect + ! + ! local workspace + ! + ! + ! max number of line segments is: + ! + ! (number of longitudes)*(max average number of crossings per line segment = 3)*ncube*2 + ! + real (kind=real_kind) , & + dimension(jmax_segments,nreconstruction), intent(out) :: weights + integer (kind=int_kind), & + dimension(jmax_segments,2), intent(out) :: weights_eul_index + + real (kind=real_kind), dimension(0:3) :: x,y + integer (kind=int_kind),dimension(0:5) :: jx_eul, jy_eul + integer (kind=int_kind) :: jsegment,i + ! + ! variables for registering crossings with Eulerian latitudes and longitudes + ! + integer (kind=int_kind) :: jcross_lat, iter + ! + ! max. crossings per side is 2*nhe + ! + real (kind=real_kind), & + dimension(jmax_segments,2) :: r_cross_lat + integer (kind=int_kind), & + dimension(jmax_segments,2) :: cross_lat_eul_index + real (kind=real_kind) , dimension(1:nvertex) :: xcell,ycell + + real (kind=real_kind) :: eps + + ldbg_global = ldbg + ldbgr = ldbg + + nc = nc_in + nhe = nhe_in + + xcell = xcell_in(1:nvertex) + ycell = ycell_in(1:nvertex) + + + ! + ! this is to avoid ill-conditioning problems + ! + eps = 1.0E-9 + + jsegment = 0 + weights = 0.0D0 + jcross_lat = 0 + ! + !********************** + ! + ! Integrate cell sides + ! + !********************** + + + IF (jx<-nhe.OR.jx>nc+1+nhe.OR.jy<-nhe.OR.jy>nc+1+nhe) THEN + WRITE(*,*) "jx,jy,-nhe,nc+1+nhe",jx,jy,-nhe,nc+1+nhe + STOP + END IF + + + call side_integral(xcell,ycell,nvertex,jsegment,jmax_segments,& + weights,weights_eul_index,nreconstruction,jx,jy,xgno,ygno,jx_min, jx_max, jy_min, jy_max,& + ngauss,gauss_weights,abscissae,& + jcross_lat,r_cross_lat,cross_lat_eul_index) + + ! + !********************** + ! + ! Do inner integrals + ! + !********************** + ! + call compute_inner_line_integrals_lat_nonconvex(r_cross_lat,cross_lat_eul_index,& + jcross_lat,jsegment,jmax_segments,xgno,jx_min, jx_max, jy_min, jy_max,& + weights,weights_eul_index,& + nreconstruction,ngauss,gauss_weights,abscissae) + ! + ! collect line-segment that reside in the same Eulerian cell + ! + if (jsegment>0) then + call collect(weights,weights_eul_index,nreconstruction,jcollect,jsegment,jmax_segments) + ! + ! DBG + ! + tmp=0.0 + do i=1,jcollect + tmp=tmp+weights(i,1) + enddo + + IF (abs(tmp)>0.01) THEN + WRITE(*,*) "sum of weights too large",tmp + !stop + END IF + IF (tmp<-1.0E-9) THEN + WRITE(*,*) "sum of weights is negative - negative area?",tmp,jx,jy + ! ldbgr=.TRUE. + !stop + !!Jinbo Xie + !!turn this off for phys grid as that of E3SM + !!Jinbo Xie + END IF + else + jcollect = 0 + end if + end subroutine compute_weights_cell + + + ! + !**************************************************************************** + ! + ! organize data and store it + ! + !**************************************************************************** + ! + subroutine collect(weights,weights_eul_index,nreconstruction,jcollect,jsegment,jmax_segments) + implicit none + integer (kind=int_kind) , intent(in) :: nreconstruction + real (kind=real_kind) , dimension(jmax_segments,nreconstruction), intent(inout) :: weights + integer (kind=int_kind), dimension(jmax_segments,2 ), intent(inout) :: weights_eul_index + integer (kind=int_kind), INTENT(OUT ) :: jcollect + integer (kind=int_kind), INTENT(IN ) :: jsegment,jmax_segments + ! + ! local workspace + ! + integer (kind=int_kind) :: imin, imax, jmin, jmax, i,j,k,h + logical :: ltmp + + real (kind=real_kind) , dimension(jmax_segments,nreconstruction) :: weights_out + integer (kind=int_kind), dimension(jmax_segments,2 ) :: weights_eul_index_out + + weights_out = 0.0D0 + weights_eul_index_out = -100 + + imin = MINVAL(weights_eul_index(1:jsegment,1)) + imax = MAXVAL(weights_eul_index(1:jsegment,1)) + jmin = MINVAL(weights_eul_index(1:jsegment,2)) + jmax = MAXVAL(weights_eul_index(1:jsegment,2)) + + ltmp = .FALSE. + + jcollect = 1 + + do j=jmin,jmax + do i=imin,imax + do k=1,jsegment + if (weights_eul_index(k,1)==i.AND.weights_eul_index(k,2)==j) then + weights_out(jcollect,1:nreconstruction) = & + weights_out(jcollect,1:nreconstruction) + weights(k,1:nreconstruction) + ltmp = .TRUE. + h = k + endif + enddo + if (ltmp) then + weights_eul_index_out(jcollect,:) = weights_eul_index(h,:) + jcollect = jcollect+1 + endif + ltmp = .FALSE. + enddo + enddo + jcollect = jcollect-1 + weights = weights_out + weights_eul_index = weights_eul_index_out + end subroutine collect + ! + !***************************************************************************************** + ! + ! + ! + !***************************************************************************************** + ! + subroutine compute_inner_line_integrals_lat(r_cross_lat,cross_lat_eul_index,& + jcross_lat,jsegment,jmax_segments,xgno,jx_min,jx_max,jy_min, jy_max,weights,weights_eul_index,& + nreconstruction,ngauss,gauss_weights,abscissae)!phl add jx_min etc. + implicit none + ! + ! for Gaussian quadrature + ! + real (kind=real_kind), dimension(ngauss), intent(in) :: gauss_weights, abscissae + ! + ! variables for registering crossings with Eulerian latitudes and longitudes + ! + integer (kind=int_kind), intent(in):: jcross_lat, jmax_segments,nreconstruction,ngauss + integer (kind=int_kind), intent(inout):: jsegment + ! + ! max. crossings per side is 2*nhe + ! + real (kind=real_kind), & + dimension(jmax_segments,2), intent(in):: r_cross_lat + integer (kind=int_kind), & + dimension(jmax_segments,2), intent(in):: cross_lat_eul_index + integer (kind=int_kind), intent(in) ::jx_min, jx_max, jy_min, jy_max + real (kind=real_kind), dimension(-nhe:nc+2+nhe), intent(in) :: xgno + real (kind=real_kind) , & + dimension(jmax_segments,nreconstruction), intent(inout) :: weights + integer (kind=int_kind), & + dimension(jmax_segments,2), intent(inout) :: weights_eul_index + real (kind=real_kind) , dimension(nreconstruction) :: weights_tmp + + integer (kind=int_kind) :: imin, imax, jmin, jmax, i,j,k, isgn, h, eul_jx, eul_jy + integer (kind=int_kind) :: idx_start_y,idx_end_y + logical :: ltmp,lcontinue + real (kind=real_kind), dimension(2) :: rstart,rend,rend_tmp + real (kind=real_kind), dimension(2) :: xseg, yseg +5 FORMAT(10e14.6) + + + if (jcross_lat>0) then + do i=MINVAL(cross_lat_eul_index(1:jcross_lat,2)),MAXVAL(cross_lat_eul_index(1:jcross_lat,2)) + ! + ! find "first" crossing with Eulerian cell i + ! + do k=1,jcross_lat + if (cross_lat_eul_index(k,2)==i) exit + enddo + do j=k+1,jcross_lat + ! + ! find "second" crossing with Eulerian cell i + ! + if (cross_lat_eul_index(j,2)==i) then + if (r_cross_lat(k,1)0) then + do i=MINVAL(cross_lat_eul_index(1:jcross_lat,2)),MAXVAL(cross_lat_eul_index(1:jcross_lat,2)) + ! WRITE(*,*) "looking at latitude ",i !xxxx + count = 1 + ! + ! find all crossings with Eulerian latitude i + ! + do k=1,jcross_lat + if (cross_lat_eul_index(k,2)==i) then + ! WRITE(*,*) "other crossings with latitude",i ," is ",k!xxxx + r_cross_lat_seg (count,:) = r_cross_lat (k,:) + cross_lat_eul_index_seg(count,:) = cross_lat_eul_index(k,:) + + IF (ldbg_global) then + WRITE(*,*) r_cross_lat_seg(count,1),r_cross_lat_seg(count,2) + WRITE(*,*) " " + END IF + count = count+1 + end if + enddo + count = count-1 + IF (ABS((count/2)-DBLE(count)/2.0)1000) THEN + WRITE(*,*) "search not converging",iter + STOP + END IF + lsame_cell_x = (x(2).GE.xgno(jx_eul).AND.x(2).LE.xgno(jx_eul+1)) + lsame_cell_y = (y(2).GE.ygno(jy_eul).AND.y(2).LE.ygno(jy_eul+1)) +! IF (ldbgr) WRITE(*,*) "lsame_cell_x,lsame_cell_y=",lsame_cell_x,lsame_cell_y + IF (lsame_cell_x.AND.lsame_cell_y) THEN + ! + !**************************** + ! + ! same cell integral + ! + !**************************** + ! +! IF (ldbgr) WRITE(*,*) "same cell integral",jx_eul,jy_eul + xseg(1) = x(1); yseg(1) = y(1); xseg(2) = x(2); yseg(2) = y(2) + jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul; + lcontinue = .FALSE. + ! + ! prepare for next side if (x(2),y(2)) is on a grid line + ! + IF (x(2).EQ.xgno(jx_eul+1).AND.x(3)>xgno(jx_eul+1)) THEN + ! + ! cross longitude jx_eul+1 + ! +! IF (ldbgr) WRITE(*,*) "cross longitude",jx_eul+1 + jx_eul=jx_eul+1 + ELSE IF (x(2).EQ.xgno(jx_eul ).AND.x(3)ygno(jy_eul+1)) THEN + ! + ! register crossing with latitude: line-segments point Northward + ! + jcross_lat = jcross_lat + 1 + jy_eul = jy_eul + 1 +! IF (ldbgr) WRITE(*,*) "cross latitude",jy_eul + cross_lat_eul_index(jcross_lat,1) = jx_eul + cross_lat_eul_index(jcross_lat,2) = jy_eul + r_cross_lat(jcross_lat,1) = x(2) + r_cross_lat(jcross_lat,2) = y(2) + ELSE IF (y(2).EQ.ygno(jy_eul ).AND.y(3)y(1) else "0" + ysgn2 = INT(SIGN(1.0D0,y(2)-y(1))) !"1" if y(2)>y(1) else "-1" + ! + !******************************************************************************* + ! + ! there is at least one crossing with latitudes but no crossing with longitudes + ! + !******************************************************************************* + ! + yeul = ygno(jy_eul+ysgn1) + IF (x(1).EQ.x(2)) THEN + ! + ! line segment is parallel to longitude (infinite slope) + ! +! IF (ldbgr) WRITE(*,*) "line segment parallel to longitude" + xcross = x(1) + ELSE + slope = (y(2)-y(1))/(x(2)-x(1)) + xcross = x_cross_eul_lat(x(1),y(1),yeul,slope) + ! + ! constrain crossing to be "physically" possible + ! + xcross = MIN(MAX(xcross,xgno(jx_eul)),xgno(jx_eul+1)) + + +! IF (ldbgr) WRITE(*,*) "cross latitude" + ! + ! debugging + ! + IF (xcross.GT.xgno(jx_eul+1).OR.xcross.LT.xgno(jx_eul)) THEN + WRITE(*,*) "xcross is out of range",jx,jy + WRITE(*,*) "xcross-xgno(jx_eul+1), xcross-xgno(jx_eul))",& + xcross-xgno(jx_eul+1), xcross-ygno(jx_eul) + STOP + END IF + END IF + xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xcross; yseg(2) = yeul + jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul; + ! + ! prepare for next iteration + ! + x(0) = x(1); y(0) = y(1); x(1) = xcross; y(1) = yeul; jy_eul = jy_eul+ysgn2 + ! + ! register crossing with latitude + ! + jcross_lat = jcross_lat+1 + cross_lat_eul_index(jcross_lat,1) = jx_eul + if (ysgn2>0) then + cross_lat_eul_index(jcross_lat,2) = jy_eul + else + cross_lat_eul_index(jcross_lat,2) = jy_eul+1 + end if + r_cross_lat(jcross_lat,1) = xcross + r_cross_lat(jcross_lat,2) = yeul + ELSE IF (lsame_cell_y) THEN +! IF (ldbgr) WRITE(*,*) "same cell y" + ! + !******************************************************************************* + ! + ! there is at least one crossing with longitudes but no crossing with latitudes + ! + !******************************************************************************* + ! + xsgn1 = (1+INT(SIGN(1.0D0,x(2)-x(1))))/2 !"1" if x(2)>x(1) else "0" + xsgn2 = INT(SIGN(1.0D0,x(2)-x(1))) !"1" if x(2)>x(1) else "-1" + xeul = xgno(jx_eul+xsgn1) +! IF (ldbgr) WRITE(*,*) " crossing longitude",jx_eul+xsgn1 + IF (ABS(x(2)-x(1))x(1) else "0" + xsgn2 = (INT(SIGN(1.0D0,x(2)-x(1)))) !"1" if x(2)>x(1) else "0" + xeul = xgno(jx_eul+xsgn1) + ysgn1 = (1+INT(SIGN(1.0D0,y(2)-y(1))))/2 !"1" if y(2)>y(1) else "0" + ysgn2 = INT(SIGN(1.0D0,y(2)-y(1))) !"1" if y(2)>y(1) else "-1" + yeul = ygno(jy_eul+ysgn1) + + slope = (y(2)-y(1))/(x(2)-x(1)) + IF (ABS(x(2)-x(1))0.AND.xcross.LE.xeul).OR.(xsgn2<0.AND.xcross.GE.xeul)) THEN + ! + ! cross latitude + ! +! IF (ldbgr) WRITE(*,*) "crossing latitude",jy_eul+ysgn1 + xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xcross; yseg(2) = yeul + jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul; + ! + ! prepare for next iteration + ! + x(0) = x(1); y(0) = y(1); x(1) = xcross; y(1) = yeul; jy_eul = jy_eul+ysgn2 + ! + ! register crossing with latitude + ! + jcross_lat = jcross_lat+1 + cross_lat_eul_index(jcross_lat,1) = jx_eul + if (ysgn2>0) then + cross_lat_eul_index(jcross_lat,2) = jy_eul + else + cross_lat_eul_index(jcross_lat,2) = jy_eul+1 + end if + r_cross_lat(jcross_lat,1) = xcross + r_cross_lat(jcross_lat,2) = yeul + ELSE + ! + ! cross longitude + ! +! IF (ldbgr) WRITE(*,*) "crossing longitude",jx_eul+xsgn1 + xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xeul; yseg(2) = ycross + jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul; + ! + ! prepare for next iteration + ! + x(0) = x(1); y(0) = y(1); x(1) = xeul; y(1) = ycross; jx_eul = jx_eul+xsgn2 + END IF + + END IF + END IF + ! + ! register line-segment (don't register line-segment if outside of panel) + ! + if (jx_eul_tmp>=jx_min.AND.jy_eul_tmp>=jy_min.AND.& + jx_eul_tmp<=jx_max-1.AND.jy_eul_tmp<=jy_max-1) then + ! jx_eul_tmp<=jx_max-1.AND.jy_eul_tmp<=jy_max-1.AND.side_count<3) then + jsegment=jsegment+1 + weights_eul_index(jsegment,1) = jx_eul_tmp + weights_eul_index(jsegment,2) = jy_eul_tmp + call get_weights_gauss(weights(jsegment,1:nreconstruction),& + xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae) + +! if (ldbg_global) then +! OPEN(unit=40, file='side_integral.dat',status='old',access='append') +! WRITE(40,*) xseg(1),yseg(1) +! WRITE(40,*) xseg(2),yseg(2) +! WRITE(40,*) " " +! CLOSE(40) +! end if + + + jdbg=jdbg+1 + + if (xseg(1).EQ.xseg(2))then + slope = bignum + else if (abs(yseg(1) -yseg(2))0) THEN + compute_slope = (y(2)-y(1))/(x(2)-x(1)) + else + compute_slope = bignum + end if + end function compute_slope + + real (kind=real_kind) function y_cross_eul_lon(x,y,xeul,slope) + implicit none + real (kind=real_kind), intent(in) :: x,y + real (kind=real_kind) , intent(in) :: xeul,slope + ! line: y=a*x+b + real (kind=real_kind) :: a,b + b = y-slope*x + y_cross_eul_lon = slope*xeul+b + end function y_cross_eul_lon + + real (kind=real_kind) function x_cross_eul_lat(x,y,yeul,slope) + implicit none + real (kind=real_kind), intent(in) :: x,y + real (kind=real_kind) , intent(in) :: yeul,slope + + if (fuzzy(ABS(slope),fuzzy_width)>0) THEN + x_cross_eul_lat = x+(yeul-y)/slope + ELSE + ! WRITE(*,*) "WARNING: slope is epsilon - ABORT" + x_cross_eul_lat = bignum + END IF + end function x_cross_eul_lat + + subroutine get_weights_exact(weights,xseg,yseg,nreconstruction) +! use cslam_analytic_mod, only: I_00, I_10, I_01, I_20, I_02, I_11 + implicit none + integer (kind=int_kind), intent(in) :: nreconstruction + real (kind=real_kind), dimension(nreconstruction), intent(out) :: weights + real (kind=real_kind), dimension(2 ), intent(in) :: xseg,yseg + ! + ! compute weights + ! + real (kind=real_kind) :: tmp,slope,b,integral,dx2,xc + integer (kind=int_kind) :: i +! weights(:) = -half*(xseg(1)*yseg(2)-xseg(2)*yseg(1)) !dummy for testing + + weights(1) = ((I_00(xseg(2),yseg(2))-I_00(xseg(1),yseg(1)))) + if (ABS(weights(1))>1.0) THEN + WRITE(*,*) "1 exact weights(jsegment)",weights(1),xseg,yseg + stop + end if + if (nreconstruction>1) then + weights(2) = ((I_10(xseg(2),yseg(2))-I_10(xseg(1),yseg(1)))) + weights(3) = ((I_01(xseg(2),yseg(2))-I_01(xseg(1),yseg(1)))) + endif + if (nreconstruction>3) then + weights(4) = ((I_20(xseg(2),yseg(2))-I_20(xseg(1),yseg(1)))) + weights(5) = ((I_02(xseg(2),yseg(2))-I_02(xseg(1),yseg(1)))) + weights(6) = ((I_11(xseg(2),yseg(2))-I_11(xseg(1),yseg(1)))) + endif + + end subroutine get_weights_exact + + + + subroutine get_weights_gauss(weights,xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae) + implicit none + integer (kind=int_kind), intent(in) :: nreconstruction,ngauss + real (kind=real_kind), dimension(nreconstruction), intent(out) :: weights + real (kind=real_kind), dimension(2 ), intent(in) :: xseg,yseg + real (kind=real_kind) :: slope + ! + ! compute weights + ! + ! + ! for Gaussian quadrature + ! + real (kind=real_kind), dimension(ngauss), intent(in) :: gauss_weights, abscissae + + ! if line-segment parallel to x or y use exact formulaes else use qudrature + ! + real (kind=real_kind) :: tmp,b,integral,dx2,xc,x,y + integer (kind=int_kind) :: i + + + + +! if (fuzzy(abs(xseg(1) -xseg(2)),fuzzy_width)==0)then + if (xseg(1).EQ.xseg(2))then + weights = 0.0D0 + else if (abs(yseg(1) -yseg(2))1) then + weights(2) = ((I_10(xseg(2),yseg(2))-I_10(xseg(1),yseg(1)))) + weights(3) = ((I_01(xseg(2),yseg(2))-I_01(xseg(1),yseg(1)))) + endif + if (nreconstruction>3) then + weights(4) = ((I_20(xseg(2),yseg(2))-I_20(xseg(1),yseg(1)))) + weights(5) = ((I_02(xseg(2),yseg(2))-I_02(xseg(1),yseg(1)))) + weights(6) = ((I_11(xseg(2),yseg(2))-I_11(xseg(1),yseg(1)))) + endif + else + + + slope = (yseg(2)-yseg(1))/(xseg(2)-xseg(1)) + b = yseg(1)-slope*xseg(1) + dx2 = 0.5D0*(xseg(2)-xseg(1)) + if (ldbgr) WRITE(*,*) "dx2 and slope in gauss weight",dx2,slope + xc = 0.5D0*(xseg(1)+xseg(2)) + integral = 0.0D0 + do i=1,ngauss + x = xc+abscissae(i)*dx2 + y = slope*x+b + integral = integral+gauss_weights(i)*F_00(x,y) + enddo + weights(1) = integral*dx2 + if (nreconstruction>1) then + integral = 0.0D0 + do i=1,ngauss + x = xc+abscissae(i)*dx2 + y = slope*x+b + integral = integral+gauss_weights(i)*F_10(x,y) + enddo + weights(2) = integral*dx2 + integral = 0.0D0 + do i=1,ngauss + x = xc+abscissae(i)*dx2 + y = slope*x+b + integral = integral+gauss_weights(i)*F_01(x,y) + enddo + weights(3) = integral*dx2 + endif + if (nreconstruction>3) then + integral = 0.0D0 + do i=1,ngauss + x = xc+abscissae(i)*dx2 + y = slope*x+b + integral = integral+gauss_weights(i)*F_20(x,y) + enddo + weights(4) = integral*dx2 + integral = 0.0D0 + do i=1,ngauss + x = xc+abscissae(i)*dx2 + y = slope*x+b + integral = integral+gauss_weights(i)*F_02(x,y) + enddo + weights(5) = integral*dx2 + integral = 0.0D0 + do i=1,ngauss + x = xc+abscissae(i)*dx2 + y = slope*x+b + integral = integral+gauss_weights(i)*F_11(x,y) + enddo + weights(6) = integral*dx2 + endif + end if + end subroutine get_weights_gauss + + real (kind=real_kind) function F_00(x_in,y_in) + implicit none + real (kind=real_kind), intent(in) :: x_in,y_in + real (kind=real_kind) :: x,y,tmp + + x = x_in + y = y_in + + F_00 =y/((1.0D0+x*x)*SQRT(1.0D0+x*x+y*y)) + end function F_00 + + real (kind=real_kind) function F_10(x_in,y_in) + implicit none + real (kind=real_kind), intent(in) :: x_in,y_in + real (kind=real_kind) :: x,y,tmp + + x = x_in + y = y_in + + F_10 =x*y/((1.0D0+x*x)*SQRT(1.0D0+x*x+y*y)) + end function F_10 + + real (kind=real_kind) function F_01(x_in,y_in) + implicit none + real (kind=real_kind), intent(in) :: x_in,y_in + real (kind=real_kind) :: x,y,tmp + + x = x_in + y = y_in + + F_01 =-1.0D0/(SQRT(1.0D0+x*x+y*y)) + end function F_01 + + real (kind=real_kind) function F_20(x_in,y_in) + implicit none + real (kind=real_kind), intent(in) :: x_in,y_in + real (kind=real_kind) :: x,y,tmp + + x = x_in + y = y_in + + F_20 =x*x*y/((1.0D0+x*x)*SQRT(1.0D0+x*x+y*y)) + end function F_20 + + real (kind=real_kind) function F_02(x_in,y_in) + implicit none + real (kind=real_kind), intent(in) :: x_in,y_in + real (kind=real_kind) :: x,y,alpha, tmp + + x = x_in + y = y_in + + alpha = ATAN(x) + tmp=y*COS(alpha) + F_02 =-y/SQRT(1.0D0+x*x+y*y)+log(tmp+sqrt(tmp*tmp+1)) + + ! + ! cos(alpha) = 1/sqrt(1+x*x) + ! + end function F_02 + + real (kind=real_kind) function F_11(x_in,y_in) + implicit none + real (kind=real_kind), intent(in) :: x_in,y_in + real (kind=real_kind) :: x,y,tmp + + x = x_in + y = y_in + + F_11 =-x/(SQRT(1.0D0+x*x+y*y)) + end function F_11 + + subroutine which_eul_cell(x,j_eul,gno) + implicit none + integer (kind=int_kind) , intent(inout) :: j_eul + real (kind=real_kind), dimension(3) , intent(in) :: x + real (kind=real_kind), dimension(-nhe:nc+2+nhe), intent(in) :: gno !phl +! real (kind=real_kind), intent(in) :: eps + + real (kind=real_kind) :: d1,d2,d3,d1p1 + logical :: lcontinue + integer :: iter + + + ! + ! this is not needed in transport code search + ! +! IF (x(1)gno(nc+2+nhe)) j_eul=nc+1+nhe +! RETURN + +! j_eul = MIN(MAX(j_eul,-nhe),nc+1+nhe) !added + + lcontinue = .TRUE. + iter = 0 + IF (ldbgr) WRITE(*,*) "from which_eul_cell",x(1),x(2),x(3) + DO WHILE (lcontinue) + iter = iter+1 + IF (x(1).GE.gno(j_eul).AND.x(1).LT.gno(j_eul+1)) THEN + lcontinue = .FALSE. + ! + ! special case when x(1) is on top of grid line + ! + IF (x(1).EQ.gno(j_eul)) THEN +! IF (ABS(x(1)-gno(j_eul))1000.OR.j_eul<-nhe.OR.j_eul>nc+2+nhe) THEN + WRITE(*,*) "search in which_eul_cell not converging!", iter,j_eul + WRITE(*,*) "input", x + WRITE(*,*) "gno", gno(nc),gno(nc+1),gno(nc+2),gno(nc+3) + STOP + END IF + END DO + END subroutine which_eul_cell + + + subroutine truncate_vertex(x,j_eul,gno) + implicit none + integer (kind=int_kind) , intent(inout) :: j_eul + real (kind=real_kind) , intent(inout) :: x + real (kind=real_kind), dimension(-nhe:nc+2+nhe), intent(in) :: gno !phl +! real (kind=real_kind), intent(in) :: eps + + logical :: lcontinue + integer :: iter + real (kind=real_kind) :: xsgn,dist,dist_new,tmp + + ! + ! this is not needed in transport code search + ! +! IF (xgno(nc+2+nhe)) j_eul=nc+1+nhe +! +! RETURN + + + lcontinue = .TRUE. + iter = 0 + dist = bignum +! j_eul = MIN(MAX(j_eul,-nhe),nc+1+nhe) !added + xsgn = INT(SIGN(1.0_dbl_kind,x-gno(j_eul))) + DO WHILE (lcontinue) + iter = iter+1 + tmp = x-gno(j_eul) + dist_new = ABS(tmp) + IF (dist_new>dist) THEN + lcontinue = .FALSE. +! ELSE IF (ABS(tmp)<1.0E-11) THEN + ELSE IF (ABS(tmp)<1.0E-9) THEN +! ELSE IF (ABS(tmp)<1.0E-4) THEN + x = gno(j_eul) + lcontinue = .FALSE. + ELSE + j_eul = j_eul+xsgn + dist = dist_new + END IF + IF (iter>10000) THEN + WRITE(*,*) "truncate vertex not converging" + STOP + END IF + END DO + END subroutine truncate_vertex + + + + +!******************************************************************************** +! +! Gauss-Legendre quadrature +! +! Tabulated values +! +!******************************************************************************** +subroutine gauss_points(n,weights,points) + implicit none + real (kind=real_kind), dimension(n), intent(out) :: weights, points + integer (kind=int_kind) , intent(in ) :: n + + select case (n) +! CASE(1) +! abscissae(1) = 0.0D0 +! weights(1) = 2.0D0 + case(2) + points(1) = -sqrt(1.0D0/3.0D0) + points(2) = sqrt(1.0D0/3.0D0) + weights(1) = 1.0D0 + weights(2) = 1.0D0 + case(3) + points(1) = -0.774596669241483377035853079956D0 + points(2) = 0.0D0 + points(3) = 0.774596669241483377035853079956D0 + weights(1) = 0.555555555555555555555555555556D0 + weights(2) = 0.888888888888888888888888888889D0 + weights(3) = 0.555555555555555555555555555556D0 + case(4) + points(1) = -0.861136311594052575223946488893D0 + points(2) = -0.339981043584856264802665659103D0 + points(3) = 0.339981043584856264802665659103D0 + points(4) = 0.861136311594052575223946488893D0 + weights(1) = 0.347854845137453857373063949222D0 + weights(2) = 0.652145154862546142626936050778D0 + weights(3) = 0.652145154862546142626936050778D0 + weights(4) = 0.347854845137453857373063949222D0 + case(5) + points(1) = -(1.0D0/3.0D0)*sqrt(5.0D0+2.0D0*sqrt(10.0D0/7.0D0)) + points(2) = -(1.0D0/3.0D0)*sqrt(5.0D0-2.0D0*sqrt(10.0D0/7.0D0)) + points(3) = 0.0D0 + points(4) = (1.0D0/3.0D0)*sqrt(5.0D0-2.0D0*sqrt(10.0D0/7.0D0)) + points(5) = (1.0D0/3.0D0)*sqrt(5.0D0+2.0D0*sqrt(10.0D0/7.0D0)) + weights(1) = (322.0D0-13.0D0*sqrt(70.0D0))/900.0D0 + weights(2) = (322.0D0+13.0D0*sqrt(70.0D0))/900.0D0 + weights(3) = 128.0D0/225.0D0 + weights(4) = (322.0D0+13.0D0*sqrt(70.0D0))/900.0D0 + weights(5) = (322.0D0-13.0D0*sqrt(70.0D0))/900.0D0 + case default + write(*,*) 'n out of range in glwp of module gll. n=',n + write(*,*) '0 0.0D0) THEN + signum = 1.0D0 + ELSEIF (x < 0.0D0) THEN + signum = -1.0D0 + ELSE + signum = 0.0D0 + ENDIF + end function + +!------------------------------------------------------------------------------ +! FUNCTION SIGNUM_FUZZY +! +! Description: +! Gives the sign of the given real number, returning zero if x is within +! a small amount from zero. +!------------------------------------------------------------------------------ + function signum_fuzzy(x) + implicit none + + real (kind=real_kind) :: signum_fuzzy + real (kind=real_kind) :: x + + IF (x > fuzzy_width) THEN + signum_fuzzy = 1.0D0 + ELSEIF (x < fuzzy_width) THEN + signum_fuzzy = -1.0D0 + ELSE + signum_fuzzy = 0.0D0 + ENDIF + end function + + function fuzzy(x,epsilon) + implicit none + + integer (kind=int_kind) :: fuzzy + real (kind=real_kind), intent(in) :: epsilon + real (kind=real_kind) :: x + + IF (ABS(x)epsilon) THEN + fuzzy = 1 + ELSE !IF (x < fuzzy_width) THEN + fuzzy = -1 + ENDIF + end function + +! +! see, e.g., http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/ +! +subroutine check_lines_cross(x1,x2,x3,x4,y1,y2,y3,y4,lcross) + implicit none + real (kind=real_kind), INTENT(IN) :: x1,x2,x3,x4,y1,y2,y3,y4 + LOGICAL, INTENT(OUT) :: lcross + ! + ! local workspace + ! + real (kind=real_kind) :: cp,tx,ty + + cp = (y4-y3)*(x2-x1)-(x4-x3)*(y2-y1) + IF (ABS(cp)-tiny.AND.tx<1.0D0+tiny.AND.& + ty>-tiny.AND.ty<1.0D0+tiny) THEN + lcross = .TRUE. + ELSE + lcross = .FALSE. +! WRITE(*,*) "not parallel but not crossing,",tx,ty + ENDIF + ENDIF +end subroutine check_lines_cross + + + REAL (KIND=dbl_kind) FUNCTION I_00(x_in,y_in) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in + REAL (KIND=dbl_kind) :: x,y + + x = x_in/aa + y = y_in/aa +! x = x_in +! y = y_in + I_00 = ATAN(x*y/SQRT(one+x*x+y*y)) + END FUNCTION I_00 + + REAL (KIND=dbl_kind) FUNCTION I_10(x_in,y_in) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in + REAL (KIND=dbl_kind) :: x,y,tmp + + x = x_in/aa + y = y_in/aa + tmp = ATAN(x) + I_10 = -ASINH(y*COS(tmp)) + ! + ! = -arcsinh(y/sqrt(1+x^2)) + ! + END FUNCTION I_10 + + REAL (KIND=dbl_kind) FUNCTION I_10_ab(alpha,beta) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta + I_10_ab = -ASINH(COS(alpha) * TAN(beta)) + END FUNCTION I_10_AB + + REAL (KIND=dbl_kind) FUNCTION I_01(x_in,y_in) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in + REAL (KIND=dbl_kind) :: x,y!,beta + + x = x_in/aa + y = y_in/aa +! beta = ATAN(y) +! I_01 = -ASINH(x*COS(beta)) + I_01 = -ASINH(x/SQRT(1+y*y)) + END FUNCTION I_01 + + REAL (KIND=dbl_kind) FUNCTION I_01_ab(alpha,beta) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta + I_01_ab = -ASINH(COS(beta) * TAN(alpha)) + END FUNCTION I_01_AB + + REAL (KIND=dbl_kind) FUNCTION I_20(x_in,y_in) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in + REAL (KIND=dbl_kind) :: x,y, tmp!,alpha,beta + + x = x_in/aa + y = y_in/aa +! alpha = aa*ATAN(x) +! beta = aa*ATAN(y) + + tmp = one+y*y + +! I_20 = y*ASINH(COS(beta)*x)+ACOS(SIN(alpha)*SIN(beta)) + I_20 = y*ASINH(x/SQRT(tmp))+ACOS(x*y/(SQRT((one+x*x)*tmp))) + END FUNCTION I_20 + + REAL (KIND=dbl_kind) FUNCTION I_20_ab(alpha,beta) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta + + I_20_ab = TAN(beta)*ASINH(COS(beta)*TAN(alpha))+ACOS(SIN(alpha)*SIN(beta)) + END FUNCTION I_20_AB + + REAL (KIND=dbl_kind) FUNCTION I_02(x_in,y_in) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in + REAL (KIND=dbl_kind) :: x,y, tmp!,alpha,beta + + x = x_in/aa + y = y_in/aa +! alpha = aa*ATAN(x) +! beta = aa*ATAN(y) + + tmp=one+x*x + + I_02 = x*ASINH(y/SQRT(tmp))+ACOS(x*y/SQRT(tmp*(1+y*y))) + END FUNCTION I_02 + + REAL (KIND=dbl_kind) FUNCTION I_02_ab(alpha,beta) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta + + I_02_ab = TAN(alpha)*ASINH(TAN(beta)*COS(alpha))+ACOS(SIN(alpha)*SIN(beta)) + END FUNCTION I_02_AB + + + REAL (KIND=dbl_kind) FUNCTION I_11(x_in,y_in) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: x_in,y_in + REAL (KIND=dbl_kind) :: x,y + + x = x_in/aa + y = y_in/aa + + I_11 = -SQRT(1+x*x+y*y) + END FUNCTION I_11 + + REAL (KIND=dbl_kind) FUNCTION I_11_ab(alpha,beta) + IMPLICIT NONE + REAL (KIND=dbl_kind), INTENT(IN) :: alpha,beta + + I_11_ab = -SQRT(one+TAN(alpha)**2+TAN(beta)**2) + END FUNCTION I_11_AB +!------------------------------------------------------------------------------ +! FUNCTION ASINH +! +! Description: +! Hyperbolic arcsin function +!------------------------------------------------------------------------------ + FUNCTION ASINH(x) + IMPLICIT NONE + + REAL (KIND=dbl_kind) :: ASINH + REAL (KIND=dbl_kind) :: x + + ASINH = LOG(x + SQRT(x * x + one)) + END FUNCTION + + + !******************************************************************************** + ! + ! Gauss-Legendre quadrature + ! + ! Tabulated values + ! + !******************************************************************************** + SUBROUTINE glwp(n,weights,abscissae) + IMPLICIT NONE + REAL (KIND=dbl_kind), DIMENSION(n), INTENT(OUT) :: weights, abscissae + INTEGER (KIND=int_kind) , INTENT(IN ) :: n + + SELECT CASE (n) + CASE(1) + abscissae(1) = 0.0 + weights(1) = 2.0 + CASE(2) + abscissae(1) = -SQRT(1.0/3.0) + abscissae(2) = SQRT(1.0/3.0) + weights(1) = 1.0 + weights(2) = 1.0 + CASE(3) + abscissae(1) = -0.774596669241483377035853079956_dbl_kind + abscissae(2) = 0.0 + abscissae(3) = 0.774596669241483377035853079956_dbl_kind + weights(1) = 0.555555555555555555555555555556_dbl_kind + weights(2) = 0.888888888888888888888888888889_dbl_kind + weights(3) = 0.555555555555555555555555555556_dbl_kind + CASE(4) + abscissae(1) = -0.861136311594052575223946488893_dbl_kind + abscissae(2) = -0.339981043584856264802665659103_dbl_kind + abscissae(3) = 0.339981043584856264802665659103_dbl_kind + abscissae(4) = 0.861136311594052575223946488893_dbl_kind + weights(1) = 0.347854845137453857373063949222_dbl_kind + weights(2) = 0.652145154862546142626936050778_dbl_kind + weights(3) = 0.652145154862546142626936050778_dbl_kind + weights(4) = 0.347854845137453857373063949222_dbl_kind + CASE(5) + abscissae(1) = -(1.0/3.0)*SQRT(5.0+2.0*SQRT(10.0/7.0)) + abscissae(2) = -(1.0/3.0)*SQRT(5.0-2.0*SQRT(10.0/7.0)) + abscissae(3) = 0.0 + abscissae(4) = (1.0/3.0)*SQRT(5.0-2.0*SQRT(10.0/7.0)) + abscissae(5) = (1.0/3.0)*SQRT(5.0+2.0*SQRT(10.0/7.0)) + weights(1) = (322.0_dbl_kind-13.0_dbl_kind*SQRT(70.0_dbl_kind))/900.0_dbl_kind + weights(2) = (322.0_dbl_kind+13.0_dbl_kind*SQRT(70.0_dbl_kind))/900.0_dbl_kind + weights(3) = 128.0_dbl_kind/225.0_dbl_kind + weights(4) = (322.0_dbl_kind+13.0_dbl_kind*SQRT(70.0_dbl_kind))/900.0_dbl_kind + weights(5) = (322.0_dbl_kind-13.0_dbl_kind*SQRT(70.0_dbl_kind))/900.0_dbl_kind + CASE DEFAULT + WRITE(*,*) 'n out of range in glwp of module gll. n=',n + WRITE(*,*) '0 shr_kind_r8 +contains +!------------------------------------------------------------------------------ +! SUBROUTINE CubedSphereABPFromRLL +! +! Description: +! Determine the (alpha,beta,panel) coordinate of a point on the sphere from +! a given regular lat lon coordinate. +! +! Parameters: +! lon - Coordinate longitude +! lat - Coordinate latitude +! alpha (OUT) - Alpha coordinate +! beta (OUT) - Beta coordinate +! ipanel (OUT) - Face panel +!------------------------------------------------------------------------------ +SUBROUTINE CubedSphereABPFromRLL(lon, lat, alpha, beta, ipanel, ldetermine_panel) + use shr_kind_mod, only: r8 => shr_kind_r8 + IMPLICIT NONE + + REAL (R8), INTENT(IN) :: lon, lat + REAL (R8), INTENT(OUT) :: alpha, beta + INTEGER :: ipanel + LOGICAL, INTENT(IN) :: ldetermine_panel + REAL (r8), PARAMETER :: pi = 3.14159265358979323846264338327 + REAL (r8), PARAMETER :: piq = 0.25*pi + REAL (r8), PARAMETER :: rotate_cube = 0.0 + + ! Local variables + REAL (R8) :: xx, yy, zz, pm + REAL (R8) :: sx, sy, sz + INTEGER :: ix, iy, iz + + ! Translate to (x,y,z) space + xx = COS(lon-rotate_cube) * COS(lat) + yy = SIN(lon-rotate_cube) * COS(lat) + zz = SIN(lat) + + pm = MAX(ABS(xx), ABS(yy), ABS(zz)) + + ! Check maximality of the x coordinate + IF (pm == ABS(xx)) THEN + IF (xx > 0) THEN; ix = 1; ELSE; ix = -1; ENDIF + ELSE + ix = 0 + ENDIF + + ! Check maximality of the y coordinate + IF (pm == ABS(yy)) THEN + IF (yy > 0) THEN; iy = 1; ELSE; iy = -1; ENDIF + ELSE + iy = 0 + ENDIF + + ! Check maximality of the z coordinate + IF (pm == ABS(zz)) THEN + IF (zz > 0) THEN; iz = 1; ELSE; iz = -1; ENDIF + ELSE + iz = 0 + ENDIF + + ! Panel assignments + IF (ldetermine_panel) THEN + IF (iz == 1) THEN + ipanel = 6; sx = yy; sy = -xx; sz = zz + + ELSEIF (iz == -1) THEN + ipanel = 5; sx = yy; sy = xx; sz = -zz + + ELSEIF ((ix == 1) .AND. (iy /= 1)) THEN + ipanel = 1; sx = yy; sy = zz; sz = xx + + ELSEIF ((ix == -1) .AND. (iy /= -1)) THEN + ipanel = 3; sx = -yy; sy = zz; sz = -xx + + ELSEIF ((iy == 1) .AND. (ix /= -1)) THEN + ipanel = 2; sx = -xx; sy = zz; sz = yy + + ELSEIF ((iy == -1) .AND. (ix /= 1)) THEN + ipanel = 4; sx = xx; sy = zz; sz = -yy + + ELSE + WRITE(*,*) 'Fatal Error: CubedSphereABPFromRLL failed' + WRITE(*,*) '(xx, yy, zz) = (', xx, ',', yy, ',', zz, ')' + WRITE(*,*) 'pm =', pm, ' (ix, iy, iz) = (', ix, ',', iy, ',', iz, ')' + STOP + ENDIF + ELSE + IF (ipanel == 6) THEN + sx = yy; sy = -xx; sz = zz + ELSEIF (ipanel == 5) THEN + sx = yy; sy = xx; sz = -zz + ELSEIF (ipanel == 1) THEN + sx = yy; sy = zz; sz = xx + ELSEIF (ipanel == 3) THEN + sx = -yy; sy = zz; sz = -xx + ELSEIF (ipanel == 2) THEN + sx = -xx; sy = zz; sz = yy + ELSEIF (ipanel == 4) THEN + sx = xx; sy = zz; sz = -yy + ELSE + WRITE(*,*) "ipanel out of range",ipanel + STOP + END IF + END IF + + ! Use panel information to calculate (alpha, beta) coords + alpha = ATAN(sx / sz) + beta = ATAN(sy / sz) + +END SUBROUTINE CubedSphereABPFromRLL + +!------------------------------------------------------------------------------ +! SUBROUTINE EquiangularAllAreas +! +! Description: +! Compute the area of all cubed sphere grid cells, storing the results in +! a two dimensional array. +! +! Parameters: +! icube - Resolution of the cubed sphere +! dA (OUT) - Output array containing the area of all cubed sphere grid cells +!------------------------------------------------------------------------------ +SUBROUTINE EquiangularAllAreas(icube, dA) + use shr_kind_mod, only: r8 => shr_kind_r8 + IMPLICIT NONE + + INTEGER, INTENT(IN) :: icube + REAL (r8), DIMENSION(icube,icube), INTENT(OUT) :: dA + + ! Local variables + INTEGER :: k, k1, k2 + REAL (r8) :: a1, a2, a3, a4 + REAL (r8), DIMENSION(icube+1,icube+1) :: ang + REAL (r8), DIMENSION(icube+1) :: gp + + REAL (r8), PARAMETER :: pi = 3.14159265358979323846264338327 + REAL (r8), PARAMETER :: piq = 0.25*pi + + + !#ifdef DBG + REAL (r8) :: dbg1 !DBG + !#endif + + ! Recall that we are using equi-angular spherical gridding + ! Compute the angle between equiangular cubed sphere projection grid lines. + DO k = 1, icube+1 + gp(k) = -piq + (pi/DBLE(2*(icube))) * DBLE(k-1) + ENDDO + + DO k2=1,icube+1 + DO k1=1,icube+1 + ang(k1,k2) =ACOS(-SIN(gp(k1)) * SIN(gp(k2))) + ENDDO + ENDDO + + DO k2=1,icube + DO k1=1,icube + a1 = ang(k1 , k2 ) + a2 = pi - ang(k1+1, k2 ) + a3 = pi - ang(k1 , k2+1) + a4 = ang(k1+1, k2+1) + ! area = r*r*(-2*pi+sum(interior angles)) + DA(k1,k2) = -2.0*pi+a1+a2+a3+a4 + ENDDO + ENDDO + + !#ifdef DBG + ! Only for debugging - test consistency + dbg1 = 0.0 !DBG + DO k2=1,icube + DO k1=1,icube + dbg1 = dbg1 + DA(k1,k2) !DBG + ENDDO + ENDDO + write(*,*) 'DAcube consistency: ',dbg1-4.0*pi/6.0 !DBG + !#endif +END SUBROUTINE EquiangularAllAreas + + +!------------------------------------------------------------------------------ +! SUBROUTINE CubedSphereRLLFromABP +! +! Description: +! Determine the lat lon coordinate of a point on a sphere given its +! (alpha,beta,panel) coordinate. +! +! Parameters: +! alpha - Alpha coordinate +! beta - Beta coordinate +! panel - Cubed sphere panel id +! lon (OUT) - Calculated longitude +! lat (OUT) - Calculated latitude +!------------------------------------------------------------------------------ +SUBROUTINE CubedSphereRLLFromABP(alpha, beta, ipanel, lon, lat) + use shr_kind_mod, only: r8 => shr_kind_r8 + IMPLICIT NONE + REAL (r8), INTENT(IN) :: alpha, beta + INTEGER , INTENT(IN) :: ipanel + REAL (r8), INTENT(OUT) :: lon, lat + ! Local variables + REAL (r8) :: xx, yy, zz, rotate_cube + REAL (r8), PARAMETER :: pi = 3.14159265358979323846264338327 + REAL (r8), PARAMETER :: piq = 0.25*pi + + rotate_cube = 0.0 + ! Convert to cartesian coordinates + CALL CubedSphereXYZFromABP(alpha, beta, ipanel, xx, yy, zz) + ! Convert back to lat lon + lat = ASIN(zz) + if (xx==0.0.and.yy==0.0) THEN + lon = 0.0 + else + lon = ATAN2(yy, xx) +rotate_cube + IF (lon<0.0) lon=lon+2.0*pi + IF (lon>2.0*pi) lon=lon-2.0*pi + end if +END SUBROUTINE CubedSphereRLLFromABP + +!------------------------------------------------------------------------------ +! SUBROUTINE CubedSphereXYZFromABP +! +! Description: +! Determine the Cartesian coordinate of a point on a sphere given its +! (alpha,beta,panel) coordinate. +! +! Parameters: +! alpha - Alpha coordinate +! beta - Beta coordinate +! panel - Cubed sphere panel id +! xx (OUT) - Calculated x coordinate +! yy (OUT) - Calculated y coordinate +! zz (OUT) - Calculated z coordinate +!------------------------------------------------------------------------------ +SUBROUTINE CubedSphereXYZFromABP(alpha, beta, ipanel, xx, yy, zz) + use shr_kind_mod, only: r8 => shr_kind_r8 + IMPLICIT NONE + + REAL (r8), INTENT(IN) :: alpha, beta + INTEGER , INTENT(IN) :: ipanel + REAL (r8), INTENT(OUT) :: xx, yy, zz + ! Local variables + REAL (r8) :: a1, b1, pm + REAL (r8) :: sx, sy, sz + + ! Convert to Cartesian coordinates + a1 = TAN(alpha) + b1 = TAN(beta) + + sz = (1.0 + a1 * a1 + b1 * b1)**(-0.5) + sx = sz * a1 + sy = sz * b1 + ! Panel assignments + IF (ipanel == 6) THEN + yy = sx; xx = -sy; zz = sz + ELSEIF (ipanel == 5) THEN + yy = sx; xx = sy; zz = -sz + ELSEIF (ipanel == 1) THEN + yy = sx; zz = sy; xx = sz + ELSEIF (ipanel == 3) THEN + yy = -sx; zz = sy; xx = -sz + ELSEIF (ipanel == 2) THEN + xx = -sx; zz = sy; yy = sz + ELSEIF (ipanel == 4) THEN + xx = sx; zz = sy; yy = -sz + ELSE + WRITE(*,*) 'Fatal Error: Panel out of range in CubedSphereXYZFromABP' + WRITE(*,*) '(alpha, beta, panel) = (', alpha, ',', beta, ',', ipanel, ')' + STOP + ENDIF +END SUBROUTINE CubedSphereXYZFromABP + + +SUBROUTINE remove_duplicates_integer(n_in,f_in,n_out,f_out) + use shr_kind_mod, only: r8 => shr_kind_r8 + integer, intent(in) :: n_in + integer,dimension(n_in), intent(in) :: f_in + integer, intent(out) :: n_out + integer,dimension(n_in), intent(out) :: f_out + ! + ! local work space + ! + integer :: k,i,j + ! + ! remove duplicates in ipanel_tmp + ! + k = 1 + f_out(1) = f_in(1) + outer: do i=2,n_in + do j=1,k + ! if (f_out(j) == f_in(i)) then + if (ABS(f_out(j)-f_in(i))<1.0E-10) then + ! Found a match so start looking again + cycle outer + end if + end do + ! No match found so add it to the output + k = k + 1 + f_out(k) = f_in(i) + end do outer + n_out = k +END SUBROUTINE remove_duplicates_integer + +SUBROUTINE remove_duplicates_latlon(n_in,lon_in,lat_in,n_out,lon_out,lat_out,tiny,ldbg) + use shr_kind_mod, only: r8 => shr_kind_r8 + integer, intent(in) :: n_in + real(r8),dimension(n_in), intent(inout) :: lon_in,lat_in + real, intent(in) :: tiny + integer, intent(out) :: n_out + real(r8),dimension(n_in), intent(out) :: lon_out,lat_out + logical :: ldbg + ! + ! local work space + ! + integer :: k,i,j + REAL (r8), PARAMETER :: pi = 3.14159265358979323846264338327 + REAL (r8), PARAMETER :: pih = 0.50*pi + ! + ! for pole points: make sure the longitudes are identical so that algorithm below works properly + ! + do i=2,n_in + if (abs(lat_in(i)-pih)