Skip to content

391 transform merge #78

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
nichannah wants to merge 4 commits into from
Closed

391 transform merge #78

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5f2a55d
Transform test. See https://github.com/NOAA-GFDL/MOM6/issues/391
nichannah Mar 10, 2017
8013dc7
Transform test. See https://github.com/NOAA-GFDL/MOM6/issues/391
nichannah Mar 10, 2017
55aa961
Fix merge. #391
Aug 5, 2017
f4a50ec
Update transform test for new coupling structure. #391
Aug 4, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
124 changes: 91 additions & 33 deletions src/SIS_debugging.F90
View file
Open in desktop

Large diffs are not rendered by default.

36 changes: 28 additions & 8 deletions src/SIS_dyn_cgrid.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ module SIS_dyn_cgrid
use MOM_hor_index, only : hor_index_type
use MOM_io, only : open_file, APPEND_FILE, ASCII_FILE, MULTIPLE, SINGLE_FILE
use SIS_hor_grid, only : SIS_hor_grid_type
use MOM_transform_test, only : do_transform_on_this_pe
use fms_io_mod, only : register_restart_field, restart_file_type
use fms_io_mod, only : restore_state, query_initialized
use mpp_domains_mod, only : domain2D
Expand Down Expand Up @@ -627,9 +628,17 @@ subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
logical :: do_hifreq_output ! If true, output occurs every iterative step.
logical :: do_trunc_its ! If true, overly large velocities in the iterations are truncated.
integer :: halo_sh_Ds ! The halo size that can be used in calculating sh_Ds.
real :: pos_sign, neg_sign
integer :: i, j, isc, iec, jsc, jec, n
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec

pos_sign = 1.0
neg_sign = -1.0
if (do_transform_on_this_pe()) then
pos_sign = -1.0
neg_sign = 1.0
endif

if (.not.associated(CS)) call SIS_error(FATAL, &
"SIS_C_dynamics: Module must be initialized before it is used.")

Expand All @@ -644,6 +653,8 @@ subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
Cor_u(:,:) = 0.0 ; Cor_v(:,:) = 0.0
fxic_d(:,:) = 0.0 ; fyic_d(:,:) = 0.0 ; fxic_t(:,:) = 0.0 ; fyic_t(:,:) = 0.0
fxic_s(:,:) = 0.0 ; fyic_s(:,:) = 0.0
PFu(:,:) = 0.0 ; PFv(:,:) = 0.0
CS%str_t(:, :) = 0.0

if (CS%SLAB_ICE) then
ui(:,:) = uo(:,:) ; vi(:,:) = vo(:,:)
Expand Down Expand Up @@ -917,10 +928,12 @@ subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
if (halo_sh_Ds < 2) call pass_var(sh_Ds, G%Domain, position=CORNER)
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,sh_Dt,sh_Dd,dy_dxT,dx_dyT,G,ui,vi)
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
sh_Dt(i,j) = (dy_dxT(i,j)*(G%IdyCu(I,j) * ui(I,j) - &
sh_Dt(i,j) = pos_sign * &
(dy_dxT(i,j)*(G%IdyCu(I,j) * ui(I,j) - &
G%IdyCu(I-1,j)*ui(I-1,j)) - &
dx_dyT(i,j)*(G%IdxCv(i,J) * vi(i,J) - &
G%IdxCv(i,J-1)*vi(i,J-1)))

sh_Dd(i,j) = (G%IareaT(i,j)*(G%dyCu(I,j) * ui(I,j) - &
G%dyCu(I-1,j)*ui(I-1,j)) + &
G%IareaT(i,j)*(G%dxCv(i,J) * vi(i,J) - &
Expand Down Expand Up @@ -1011,12 +1024,14 @@ subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
! Save the current values of u for later use in updating v.
u_tmp(I,j) = ui(I,j)

Cor = ((azon(I,j) * vi(i+1,J) + czon(I,j) * vi(i,J-1)) + &
(bzon(I,j) * vi(i,J) + dzon(I,j) * vi(i+1,J-1))) ! - Cor_ref_u(I,j)
Cor = pos_sign * &
((azon(I,j) * vi(i+1,J) + czon(I,j) * vi(i,J-1)) + &
(bzon(I,j) * vi(i,J) + dzon(I,j) * vi(i+1,J-1))) ! - Cor_ref_u(I,j)
! Evaluate 1/m x.Div(m strain). This expressions include all metric terms
! for an orthogonal grid. The str_d term integrates out to no curl, while
! str_s & str_t terms impose no divergence and do not act on solid body rotation.
fxic_now = G%IdxCu(I,j) * (CS%str_d(i+1,j) - CS%str_d(i,j)) + &
fxic_now = pos_sign * &
G%IdxCu(I,j) * (CS%str_d(i+1,j) - CS%str_d(i,j)) + &
(G%IdyCu(I,j)*(dy2T(i+1,j)*CS%str_t(i+1,j) - &
dy2T(i,j) *CS%str_t(i,j)) + &
G%IdxCu(I,j)*(dx2B(I,J) *CS%str_s(I,J) - &
Expand Down Expand Up @@ -1093,13 +1108,16 @@ subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
!$OMP private(Cor,fyic_now,u2_at_v,vio_init,drag_v, &
!$OMP m_vio_explicit,b_vel0,vio_pred,vio_C)
do J=jsc-1,jec ; do i=isc,iec
Cor = -1.0*((amer(I-1,j) * u_tmp(I-1,j) + cmer(I,j+1) * u_tmp(I,j+1)) + &
(bmer(I,j) * u_tmp(I,j) + dmer(I-1,j+1) * u_tmp(I-1,j+1)))

Cor = neg_sign * &
((amer(I-1,j) * u_tmp(I-1,j) + cmer(I,j+1) * u_tmp(I,j+1)) + &
(bmer(I,j) * u_tmp(I,j) + dmer(I-1,j+1) * u_tmp(I-1,j+1)))
! Evaluate 1/m y.Div(m strain). This expressions include all metric terms
! for an orthogonal grid. The str_d term integrates out to no curl, while
! str_s & str_t terms impose no divergence and do not act on solid body rotation.
fyic_now = G%IdyCv(i,J) * (CS%str_d(i,j+1)-CS%str_d(i,j)) + &
(-G%IdxCv(i,J)*(dx2T(i,j+1)*CS%str_t(i,j+1) - &
fyic_now = neg_sign * &
G%IdyCv(i,J) * (CS%str_d(i,j+1)-CS%str_d(i,j)) + &
(G%IdxCv(i,J)*(dx2T(i,j+1)*CS%str_t(i,j+1) - &
dx2T(i,j) *CS%str_t(i,j)) + &
G%IdyCv(i,J)*(dy2B(I,J) *CS%str_s(I,J) - &
dy2B(I-1,J)*CS%str_s(I-1,J)) )*G%IareaCv(i,J)
Expand Down Expand Up @@ -1207,6 +1225,8 @@ subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
call uvchksum("f[xy]ic in SIS_C_dynamics", fxic, fyic, G)
call uvchksum("f[xy]oc in SIS_C_dynamics", fxoc, fyoc, G)
call uvchksum("Cor_[uv] in SIS_C_dynamics", Cor_u, Cor_v, G)
call uvchksum("f[xy]oc in SIS_C_dynamics", fxoc, fyoc, G)
call uvchksum("f[xy]ic in SIS_C_dynamics", fxic, fyic, G)
call uvchksum("[uv]i in SIS_C_dynamics", ui, vi, G)
endif

Expand Down
6 changes: 5 additions & 1 deletion src/SIS_dyn_trans.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -315,6 +315,9 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; NkIce = IG%NkIce
Idt_slow = 0.0 ; if (dt_slow > 0.0) Idt_slow = 1.0/dt_slow

WindStr_x_Cu(:, :) = 0.0
WindStr_y_Cv(:, :) = 0.0

if (CS%specified_ice) then
ndyn_steps = 0 ; dt_slow_dyn = 0.0
!$OMP parallel do default(none) shared(isd,ied,jsd,jed,WindStr_x_A,WindStr_y_A, &
Expand Down Expand Up @@ -464,8 +467,9 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
call hchksum(OSS%sea_lev, "sea_lev before SIS_C_dynamics", G%HI, haloshift=1)
call hchksum(ice_cover, "ice_cover before SIS_C_dynamics", G%HI, haloshift=1)
call uvchksum("[uv]_ocn before SIS_C_dynamics", OSS%u_ocn_C, OSS%v_ocn_C, G, halos=1)
call uvchksum("WindStr_[xy] before SIS_C_dynamics", WindStr_x_Cu, WindStr_y_Cv, G, halos=1)
call hchksum_pair("FIA%WindStr_[xy] before SIS_C_dynamics", FIA%WindStr_x, FIA%WindStr_y, G, halos=1)
call hchksum_pair("WindStr_[xy]_A before SIS_C_dynamics", WindStr_x_A, WindStr_y_A, G, halos=1)
call uvchksum("WindStr_[xy]_C[uv] before SIS_C_dynamics", WindStr_x_Cu, WindStr_y_Cv, G, halos=1)
endif

call mpp_clock_begin(iceClocka)
Expand Down
26 changes: 20 additions & 6 deletions src/SIS_fast_thermo.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,13 @@
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
module SIS_fast_thermo

use mpp_mod, only: mpp_pe
use SIS_diag_mediator, only : SIS_diag_ctrl
! ! use SIS_diag_mediator, only : enable_SIS_averaging, disable_SIS_averaging
! ! use SIS_diag_mediator, only : query_SIS_averaging_enabled, post_SIS_data
! ! use SIS_diag_mediator, only : register_diag_field=>register_SIS_diag_field

use SIS_debugging, only : hchksum
use SIS_debugging, only : hchksum, hchksum_pair
use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg
use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
Expand All @@ -46,11 +47,13 @@ module SIS_fast_thermo
use MOM_time_manager, only : time_type, time_type_to_real
use MOM_time_manager, only : set_date, set_time, operator(+), operator(-)
use MOM_time_manager, only : operator(>), operator(*), operator(/), operator(/=)
use MOM_transform_test, only : do_transform_on_this_pe

use coupler_types_mod, only : coupler_3d_bc_type
use SIS_optics, only : ice_optics_SIS2, bright_ice_temp, SIS_optics_CS
use SIS_optics, only : VIS_DIR, VIS_DIF, NIR_DIR, NIR_DIF
use SIS_types, only : ice_state_type, IST_chksum, IST_bounds_check

use SIS_types, only : fast_ice_avg_type, ice_rad_type, simple_OSS_type, total_sfc_flux_type
use SIS_types, only : FIA_chksum

Expand Down Expand Up @@ -250,7 +253,7 @@ subroutine avg_top_quantities(FIA, Rad, IST, G, IG)
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(in) :: IG

real :: u, v, divid, sign
real :: u, v, divid, sign, pos_sign
real :: I_avc ! The inverse of the number of contributions.
real :: I_wts ! 1.0 / ice_cover or 0 if ice_cover is 0, nondim.
integer :: i, j, k, m, n, b, nb, isc, iec, jsc, jec, ncat
Expand All @@ -269,15 +272,16 @@ subroutine avg_top_quantities(FIA, Rad, IST, G, IG)
! Rotate the stress from lat/lon to ocean coordinates and possibly change the
! sign to positive for downward fluxes of positive momentum.
sign = 1.0 ; if (FIA%atmos_winds) sign = -1.0
pos_sign = 1.0 ; if (do_transform_on_this_pe()) pos_sign = -1.0
I_avc = 1.0/real(FIA%avg_count)

!$OMP parallel do default(shared) private(u,v)
do j=jsc,jec
do k=0,ncat ; do i=isc,iec
u = FIA%flux_u_top(i,j,k) * (sign*I_avc)
v = FIA%flux_v_top(i,j,k) * (sign*I_avc)
FIA%flux_u_top(i,j,k) = u*G%cos_rot(i,j)-v*G%sin_rot(i,j) ! rotate stress from lat/lon
FIA%flux_v_top(i,j,k) = v*G%cos_rot(i,j)+u*G%sin_rot(i,j) ! to ocean coordinates
FIA%flux_u_top(i,j,k) = u*G%cos_rot(i,j)-v*G%sin_rot(i,j)*pos_sign ! rotate stress from lat/lon
FIA%flux_v_top(i,j,k) = v*G%cos_rot(i,j)+u*G%sin_rot(i,j)*pos_sign ! to ocean coordinates
FIA%flux_sh_top(i,j,k) = FIA%flux_sh_top(i,j,k) * I_avc
FIA%evap_top(i,j,k) = FIA%evap_top(i,j,k) * I_avc
do b=1,nb ; FIA%flux_sw_top(i,j,k,b) = FIA%flux_sw_top(i,j,k,b) * I_avc ; enddo
Expand Down Expand Up @@ -686,6 +690,8 @@ subroutine do_update_ice_model_fast(Atmos_boundary, IST, sOSS, Rad, FIA, &
endif ; enddo ; enddo ; enddo
endif

flux_u(:, :, :) = 0.0

if (CS%debug_fast) &
call IST_chksum("Start do_update_ice_model_fast", IST, G, IG)

Expand Down Expand Up @@ -723,8 +729,7 @@ subroutine do_update_ice_model_fast(Atmos_boundary, IST, sOSS, Rad, FIA, &
enddo

if (CS%debug_fast) then
call hchksum(flux_u(:,:,1:), "Mid do_fast flux_u", G%HI)
call hchksum(flux_v(:,:,1:), "Mid do_fast flux_v", G%HI)
call hchksum_pair("Mid do_fast flux_[uv]", flux_u(:,:,1:), flux_v(:,:,1:), G%HI)
call hchksum(flux_sh(:,:,1:), "Mid do_fast flux_sh", G%HI)
call hchksum(evap(:,:,1:), "Mid do_fast evap", G%HI)
call hchksum(flux_lw(:,:,1:), "Mid do_fast flux_lw", G%HI)
Expand All @@ -734,6 +739,15 @@ subroutine do_update_ice_model_fast(Atmos_boundary, IST, sOSS, Rad, FIA, &
enddo
call hchksum(lprec(:,:,1:), "Mid do_fast lprec", G%HI)
call hchksum(fprec(:,:,1:), "Mid do_fast fprec", G%HI)
call hchksum(dshdt(:,:,2:), "Mid do_fast dshdt", G%HI)
call hchksum(devapdt(:,:,1:), "Mid do_fast devapdt", G%HI)
call hchksum(dlwdt(:,:,1:), "Mid do_fast dlwdt", G%HI)
call hchksum(flux_sw(:,:,1:,nir_dir), "Mid do_fast flux_sw_nir_dir", G%HI)
call hchksum(flux_sw(:,:,1:,nir_dif), "Mid do_fast flux_sw_nir_dif", G%HI)
call hchksum(flux_sw(:,:,1:,vis_dir), "Mid do_fast flux_sw_vis_dir", G%HI)
call hchksum(flux_sw(:,:,1:,vis_dif), "Mid do_fast flux_sw_vis_dif", G%HI)
call hchksum(lprec(:,:,1:), "Mid do_fast lprec", G%HI)
call hchksum(fprec(:,:,1:), "Mid do_fast fprec", G%HI)
call hchksum(dshdt(:,:,1:), "Mid do_fast dshdt", G%HI)
call hchksum(devapdt(:,:,1:), "Mid do_fast devapdt", G%HI)
call hchksum(dlwdt(:,:,1:), "Mid do_fast dlwdt", G%HI)
Expand Down
15 changes: 0 additions & 15 deletions src/SIS_fixed_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,15 +71,6 @@ subroutine SIS_initialize_fixed(G, PF, write_geom, output_dir)
! Initialize the various masks and any masked metrics.
call initialize_masks(G, PF)

if (debug) then
call hchksum(G%bathyT, 'SIS_initialize_fixed: depth ', G%HI, &
haloshift=min(1, G%ied-G%iec, G%jed-G%jec))
call hchksum(G%mask2dT, 'SIS_initialize_fixed: mask2dT ', G%HI)
call uvchksum('SIS_initialize_fixed: mask2dC[uv] ', &
G%mask2dCu, G%mask2dCv, G)
call Bchksum(G%mask2dBu, 'SIS_initialize_fixed: mask2dBu ', G%HI)
endif

! Modulate geometric scales according to geography.
call get_param(PF, mdl, "CHANNEL_CONFIG", config, &
"A parameter that determines which set of channels are \n"//&
Expand Down Expand Up @@ -108,12 +99,6 @@ subroutine SIS_initialize_fixed(G, PF, write_geom, output_dir)
call MOM_initialize_rotation(G%CoriolisBu, G, PF)
! Calculate the components of grad f (beta)
call MOM_calculate_grad_Coriolis(G%dF_dx, G%dF_dy, G)
if (debug) then
call Bchksum(G%CoriolisBu, "SIS_initialize_fixed: f ", G%HI)
call hchksum(G%dF_dx, "SIS_initialize_fixed: dF_dx ", G%HI)
call hchksum(G%dF_dy, "SIS_initialize_fixed: dF_dy ", G%HI)
endif

call initialize_grid_rotation_angle(G, PF)

! Write out all of the grid data used by this run.
Expand Down
78 changes: 77 additions & 1 deletion src/SIS_hor_grid.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,22 +7,28 @@ module SIS_hor_grid
use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_data_domain
use mpp_domains_mod, only : mpp_get_global_domain

use MOM_transform_test, only : do_transform_on_this_pe
use MOM_hor_index, only : hor_index_type, hor_index_init
use MOM_domains, only : MOM_domain_type, get_domain_extent, compute_block_extent
use MOM_domains, only : MOM_domains_init, clone_MOM_domain
use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
use MOM_transform_test, only : transform

implicit none ; private

#include <SIS2_memory.h>

public :: set_hor_grid, SIS_hor_grid_end, set_derived_SIS_metrics, set_first_direction, isPointInCell
public :: transform_sis_hor_grid

type, public :: SIS_hor_grid_type
type(MOM_domain_type), pointer :: Domain => NULL()
type(MOM_domain_type), pointer :: Domain_aux => NULL() ! A non-symmetric auxiliary domain type.
type(hor_index_type) :: HI

type(SIS_hor_grid_type), pointer :: self_untrans

integer :: isc, iec, jsc, jec ! The range of the computational domain indices
integer :: isd, ied, jsd, jed ! and data domain indices at tracer cell centers.
integer :: isg, ieg, jsg, jeg ! The range of the global domain tracer cell indices.
Expand Down Expand Up @@ -281,6 +287,71 @@ subroutine set_hor_grid(G, param_file, HI, global_indexing)

end subroutine set_hor_grid

subroutine transform_sis_hor_grid(G, G_trans)
type(SIS_hor_grid_type), intent(in) :: G !< The horizontal grid type
type(SIS_hor_grid_type), intent(out) :: G_trans !< The horizontal grid type

call transform(G%dyT, G_trans%dxT)
call transform(G%dxT, G_trans%dyT)
call transform(G%IdyT, G_trans%IdxT)
call transform(G%IdxT, G_trans%IdyT)

call transform(G%dyBu, G_trans%dxBu)
call transform(G%dxBu, G_trans%dyBu)
call transform(G%IdyBu, G_trans%IdxBu)
call transform(G%IdxBu, G_trans%IdyBu)

call transform(G%dyCu, G_trans%dxCv)
call transform(G%dxCu, G_trans%dyCv)
call transform(G%IdyCu, G_trans%IdxCv)
call transform(G%IdxCu, G_trans%IdyCv)

call transform(G%dyCv, G_trans%dxCu)
call transform(G%dxCv, G_trans%dyCu)
call transform(G%IdyCv, G_trans%IdxCu)
call transform(G%IdxCv, G_trans%IdyCu)

call transform(G%dy_Cu, G_trans%dx_Cv)
call transform(G%dy_Cu_obc, G_trans%dx_Cv_obc)

call transform(G%dx_Cv, G_trans%dy_Cu)
call transform(G%dx_Cv_obc, G_trans%dy_Cu_obc)

call transform(G%mask2dT, G_trans%mask2dT)
call transform(G%mask2dBu, G_trans%mask2dBu)

call transform(G%mask2dCu, G_trans%mask2dCv)
call transform(G%mask2dCv, G_trans%mask2dCu)

call transform(G%geoLatT, G_trans%geoLatT)
call transform(G%geoLonT, G_trans%geoLonT)
call transform(G%geoLatBu, G_trans%geoLatBu)
call transform(G%geoLonBu, G_trans%geoLonBu)

call transform(G%geoLatCu, G_trans%geoLatCv)
call transform(G%geoLonCu, G_trans%geoLonCv)
call transform(G%geoLatCv, G_trans%geoLatCu)
call transform(G%geoLonCv, G_trans%geoLonCu)

call transform(G%areaT, G_trans%areaT)
call transform(G%IareaT, G_trans%IareaT)
call transform(G%areaBu, G_trans%areaBu)
call transform(G%IareaBu, G_trans%IareaBu)
call transform(G%areaCu, G_trans%areaCv)
call transform(G%IareaCu, G_trans%IareaCv)
call transform(G%areaCv, G_trans%areaCu)
call transform(G%IareaCv, G_trans%IareaCu)

call transform(G%sin_rot, G_trans%sin_rot)
call transform(G%cos_rot, G_trans%cos_rot)
call transform(G%bathyT, G_trans%bathyT)

call transform(G%CoriolisBu, G_trans%CoriolisBu)
call transform(G%dF_dx, G_trans%dF_dy)
call transform(G%dF_dy, G_trans%dF_dx)

end subroutine


!> set_derived_SIS_metrics calculates metric terms that are derived from other metrics.
subroutine set_derived_SIS_metrics(G)
Expand Down Expand Up @@ -371,7 +442,12 @@ subroutine set_first_direction(G, y_first)
type(SIS_hor_grid_type), intent(inout) :: G
integer, intent(in) :: y_first

G%first_direction = y_first
if (do_transform_on_this_pe()) then
G%first_direction = y_first + 1
else
G%first_direction = y_first
endif

end subroutine set_first_direction

!---------------------------------------------------------------------
Expand Down
Loading