diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index 5ff967b2c4..5890357061 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -5625,6 +5625,12 @@ Give the user the ability to input LTR high-latitude electric potential. Default: FALSE + +For 2-way coupling WACCM-X with MAGE magnetosphere +Default: FALSE + + List of full pathnames of LTR electic potential inputs for both hemispheres. diff --git a/src/ionosphere/waccmx/dpie_coupling.F90 b/src/ionosphere/waccmx/dpie_coupling.F90 index fad558296b..57542a7bf0 100644 --- a/src/ionosphere/waccmx/dpie_coupling.F90 +++ b/src/ionosphere/waccmx/dpie_coupling.F90 @@ -13,6 +13,7 @@ module dpie_coupling use perf_mod, only: t_startf, t_stopf use amie_module, only: getamie use ltr_module, only: getltr + use mage_module, only: mage_advance use edyn_solve, only: phihm use edyn_params, only: dtr, rtd use aurora_params, only: prescribed_period ! turns on overwrite of energy fields in aurora phys @@ -119,7 +120,7 @@ subroutine d_pie_init( edyn_active_in, oplus_xport_in, oplus_nsplit_in, crit_col end subroutine d_pie_init !----------------------------------------------------------------------- - subroutine d_pie_epotent( highlat_potential_model, crit_out, cols, cole, efx_phys, kev_phys, amie_in, ltr_in ) + subroutine d_pie_epotent( highlat_potential_model, crit_out, cols, cole, efx_phys, kev_phys, amie_in, ltr_in, mage_in ) use edyn_solve, only: pfrac ! NH fraction of potential (nmlonp1,nmlat0) use time_manager, only: get_curr_date use heelis, only: heelis_model @@ -142,6 +143,7 @@ subroutine d_pie_epotent( highlat_potential_model, crit_out, cols, cole, efx_phy integer, optional, intent(in) :: cols, cole logical, optional,intent(in) :: amie_in logical, optional,intent(in) :: ltr_in + logical, optional,intent(in) :: mage_in ! Prescribed energy flux real(r8), optional, intent(out) :: efx_phys(:) @@ -151,7 +153,7 @@ subroutine d_pie_epotent( highlat_potential_model, crit_out, cols, cole, efx_phy ! ! local vars ! - logical :: amie_inputs, ltr_inputs + logical :: amie_inputs, ltr_inputs, mage_inputs real(r8) :: secs ! time of day in seconds integer :: iyear,imo,iday,tod ! tod is time-of-day in seconds @@ -202,10 +204,12 @@ subroutine d_pie_epotent( highlat_potential_model, crit_out, cols, cole, efx_phy amie_inputs=.false. ltr_inputs=.false. + mage_inputs=.false. if (present(amie_in)) amie_inputs=amie_in if (present(ltr_in)) ltr_inputs= ltr_in + if (present(ltr_in)) mage_inputs= mage_in - prescribed_inputs: if (amie_inputs .or. ltr_inputs) then + prescribed_inputs: if (amie_inputs .or. ltr_inputs .or. mage_inputs) then if (.not. (present(kev_phys).and.present(efx_phys)) ) then call endrun('d_pie_epotent: kev_phys and efx_phys must be present') @@ -224,6 +228,14 @@ subroutine d_pie_epotent( highlat_potential_model, crit_out, cols, cole, efx_phy write(iulog,"('After Calling getamie >>> iamie = ', i2)") iamie end if prescribed_period = iamie == 1 + else if (mage_inputs) then + if (masterproc) then + write(iulog,*) 'Calling mage_advance >>> ' + end if + + call mage_advance( iyear, imo, iday, tod, & + prescr_phihm, prescr_efxm, prescr_kevm ) + prescribed_period = .TRUE. else if (masterproc) then write(iulog,*) 'Calling getltr >>> ' diff --git a/src/ionosphere/waccmx/edynamo.F90 b/src/ionosphere/waccmx/edynamo.F90 index be99c7dff9..df518977f4 100644 --- a/src/ionosphere/waccmx/edynamo.F90 +++ b/src/ionosphere/waccmx/edynamo.F90 @@ -32,8 +32,10 @@ module edynamo real(r8), allocatable, dimension(:,:) :: & zigm11, & ! sigma11*cos(theta0) zigmc, & ! sigmac - zigm1, & ! for Hall conductance diagnostic + zigm1, & ! for Hall conductance diagnostic (folded) + azigm1, & ! for Hall conductance diagnostic (not folded) zigm2, & ! sigma2 + azigm2, & ! for Ped conductance diagnostic (not folded) zigm22, & ! sigma22/cos(theta0) rim1,rim2, & ! see description in comment below rhs, & ! right-hand side of PDE @@ -105,7 +107,7 @@ module edynamo logical, public :: debug_hist = .false. public :: alloc_edyn, ed1, ed2, ed1_glb, ed2_glb - public :: zigm11, zigmc, zigm2, zigm22, rim1, rim2 + public :: zigm11, zigmc, zigm2, zigm22, rim1, rim2, azigm1, azigm2 public :: dynamo contains @@ -368,9 +370,15 @@ subroutine alloc_edyn allocate(zigm1(mlon00:mlon11,mlat00:mlat11) ,stat=istat) if (istat /= 0) call endrun('alloc_edyn: zigm1') zigm1 = finit + allocate(azigm1(mlon00:mlon11,mlat00:mlat11) ,stat=istat) + if (istat /= 0) call endrun('alloc_edyn: azigm1') + azigm1 = finit allocate(zigm2(mlon00:mlon11,mlat00:mlat11) ,stat=istat) if (istat /= 0) call endrun('alloc_edyn: zigm2') zigm2 = finit + allocate(azigm2(mlon00:mlon11,mlat00:mlat11) ,stat=istat) + if (istat /= 0) call endrun('alloc_edyn: azigm2') + azigm2 = finit allocate(zigm22(mlon00:mlon11,mlat00:mlat11),stat=istat) if (istat /= 0) call endrun('alloc_edyn: zigm22') zigm22 = finit @@ -528,7 +536,9 @@ subroutine fieldline_integrals( ped_mag, hal_mag, adotv1_mag, adotv2_mag, & zigm11 = finit zigm22 = finit zigm1 = finit + azigm1 = finit zigm2 = finit + azigm2 = finit zigmc = finit rim1 = finit rim2 = finit @@ -929,6 +939,10 @@ subroutine complete_integrals fmsub(:,:,6) = rim2 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,7) = zigm1 (mlon0:mlon1,mlat0:mlat1) +! Store zigm1 for coupling before folding hemispheres + azigm1(mlon0:mlon1,mlat0:mlat1) = zigm1(mlon0:mlon1,mlat0:mlat1) + azigm2(mlon0:mlon1,mlat0:mlat1) = zigm2(mlon0:mlon1,mlat0:mlat1) + call mp_mag_foldhem(fmsub,mlon0,mlon1,mlat0,mlat1,nf2d) call mp_mag_periodic_f2d(fmsub,mlon0,mlon1,mlat0,mlat1,nf2d) diff --git a/src/ionosphere/waccmx/ionosphere_interface.F90 b/src/ionosphere/waccmx/ionosphere_interface.F90 index 59e80260f0..46ec0bd3ff 100644 --- a/src/ionosphere/waccmx/ionosphere_interface.F90 +++ b/src/ionosphere/waccmx/ionosphere_interface.F90 @@ -46,6 +46,10 @@ module ionosphere_interface type(var_desc_t) :: Optm1_vdesc logical :: opmmrtm1_initialized + type(var_desc_t) :: azigm1_vdesc + type(var_desc_t) :: azigm2_vdesc + + integer :: index_ped, index_hall, index_te, index_ti integer :: index_ui, index_vi, index_wi @@ -84,10 +88,10 @@ module ionosphere_interface character(len=cl) :: amiesh_files(max_num_files) = 'NONE' character(len=cl) :: ltr_files(max_num_files) = 'NONE' - character(len=16) :: ionos_epotential_model = 'none' logical :: ionos_epotential_amie = .false. logical :: ionos_epotential_ltr = .false. + logical :: ionos_epotential_mage = .false. integer :: indxefx=-1, indxkev=-1 integer :: oplus_nlon, oplus_nlat ! Oplus grid @@ -125,6 +129,7 @@ subroutine ionosphere_readnl( nlfile ) namelist /ionosphere_nl/ oplus_adiff_limiter, oplus_shapiro_const, oplus_enforce_floor, oplus_ring_polar_filter namelist /ionosphere_nl/ ionos_epotential_model, ionos_epotential_amie, ionos_epotential_ltr, wei05_coefs_file namelist /ionosphere_nl/ amienh_files, amiesh_files, wei05_coefs_file, ltr_files + namelist /ionosphere_nl/ ionos_epotential_mage namelist /ionosphere_nl/ epot_crit_colats namelist /ionosphere_nl/ ionos_npes namelist /ionosphere_nl/ oplus_grid, edyn_grid @@ -156,6 +161,7 @@ subroutine ionosphere_readnl( nlfile ) call mpi_bcast(ionos_epotential_model, len(ionos_epotential_model), mpi_character, masterprocid, mpicom, ierr) call mpi_bcast(ionos_epotential_amie,1, mpi_logical, masterprocid, mpicom, ierr) call mpi_bcast(ionos_epotential_ltr,1, mpi_logical, masterprocid, mpicom, ierr) + call mpi_bcast(ionos_epotential_mage,1, mpi_logical, masterprocid, mpicom, ierr) call mpi_bcast(wei05_coefs_file, len(wei05_coefs_file), mpi_character, masterprocid, mpicom, ierr) call mpi_bcast(amienh_files, max_num_files*len(amienh_files(1)), mpi_character, masterprocid, mpicom, ierr) call mpi_bcast(amiesh_files, max_num_files*len(amiesh_files(1)), mpi_character, masterprocid, mpicom, ierr) @@ -197,6 +203,7 @@ subroutine ionosphere_readnl( nlfile ) write(iulog,*) 'ionosphere_readnl: ionos_epotential_model = ', trim(ionos_epotential_model) write(iulog,*) 'ionosphere_readnl: ionos_epotential_amie = ', ionos_epotential_amie write(iulog,*) 'ionosphere_readnl: ionos_epotential_ltr = ', ionos_epotential_ltr + write(iulog,*) 'ionosphere_readnl: ionos_epotential_mage = ', ionos_epotential_mage write(iulog,'(a,2(g12.4))') & 'ionosphere_readnl: epot_crit_colats = ', epot_crit_colats write(iulog,'(a,i0)') 'ionosphere_readnl: ionos_npes = ',ionos_npes @@ -234,6 +241,7 @@ subroutine ionosphere_init() use ref_pres, only: pref_edge ! target ailev(pverp) interface levels use amie_module, only: init_amie use ltr_module, only: init_ltr + use mage_module, only: mage_init use wei05sc, only: weimer05_init use phys_control, only: phys_getopts @@ -243,7 +251,7 @@ subroutine ionosphere_init() call phys_getopts(state_debug_checks_out=state_debug_checks) - if ( ionos_epotential_amie .or. ionos_epotential_ltr) then + if ( ionos_epotential_amie .or. ionos_epotential_ltr .or. ionos_epotential_mage) then call pbuf_add_field('AUREFX', 'global', dtype_r8, (/pcols/), indxefx) ! Prescribed Energy flux call pbuf_add_field('AURKEV', 'global', dtype_r8, (/pcols/), indxkev) ! Prescribed Mean energy end if @@ -380,6 +388,11 @@ subroutine ionosphere_init() call addfld ('ltr_efx_phys', horiz_only, 'I', 'mW/m2', 'LTR energy flux') call addfld ('ltr_kev_phys', horiz_only, 'I', 'keV', 'LTR mean energy') end if + if ( ionos_epotential_mage ) then + call mage_init() + call addfld ('mage_efx_phys', horiz_only, 'I', 'mW/m2', 'MAGE energy flux') + call addfld ('mage_kev_phys', horiz_only, 'I', 'keV', 'MAGE mean energy') + end if if ( trim(ionos_epotential_model) == 'weimer' ) then call weimer05_init(wei05_coefs_file) end if @@ -390,6 +403,9 @@ subroutine ionosphere_init() call addfld ('Z3GMI', (/ 'lev' /), 'I', 'm', & 'Geometric height (Interfaces)', gridname='physgrid') + call addfld ('Early_AZIGM2', horiz_only, 'I', ' ','EDYN_ZIGM2' ,gridname='gmag_grid') + call addfld ('Early_AZIGM1', horiz_only, 'I', ' ','EDYN_ZIGM2' ,gridname='gmag_grid') + end subroutine ionosphere_init !---------------------------------------------------------------------------- @@ -397,6 +413,9 @@ end subroutine ionosphere_init subroutine ionosphere_run1(pbuf2d) use physics_buffer, only: physics_buffer_desc use cam_history, only: outfld, write_inithist + use edynamo, only: azigm1, azigm2 + use savefield_waccm, only: savefld_waccm + use edyn_mpi, only: mlon0,omlon1, mlat0,mlat1 ! args type(physics_buffer_desc), pointer :: pbuf2d(:,:) @@ -412,6 +431,9 @@ subroutine ionosphere_run1(pbuf2d) real(r8), pointer :: prescr_efx(:) ! prescribed energy flux real(r8), pointer :: prescr_kev(:) ! prescribed characteristic mean energy + call savefld_waccm(azigm2(mlon0:omlon1,mlat0:mlat1),'Early_AZIGM2',1, mlon0,omlon1,mlat0,mlat1) + call savefld_waccm(azigm1(mlon0:omlon1,mlat0:mlat1),'Early_AZIGM1',1, mlon0,omlon1,mlat0,mlat1) + if( write_inithist() .and. ionos_xport_active ) then do lchnk = begchunk, endchunk call outfld ('OpTM1&IC', opmmrtm1_phys(:,:,lchnk), pcols, lchnk) @@ -420,7 +442,7 @@ subroutine ionosphere_run1(pbuf2d) nullify(prescr_efx) nullify(prescr_kev) - prescribed_epot: if ( ionos_epotential_amie .or. ionos_epotential_ltr ) then + prescribed_epot: if ( ionos_epotential_amie .or. ionos_epotential_ltr .or. ionos_epotential_mage) then blksize = 0 do lchnk = begchunk, endchunk blksize = blksize + get_ncols_p(lchnk) @@ -432,7 +454,7 @@ subroutine ionosphere_run1(pbuf2d) ! data assimilated potential call d_pie_epotent(ionos_epotential_model, epot_crit_colats, & cols=1, cole=blksize, efx_phys=prescr_efx, kev_phys=prescr_kev, & - amie_in=ionos_epotential_amie, ltr_in=ionos_epotential_ltr ) + amie_in=ionos_epotential_amie, ltr_in=ionos_epotential_ltr, mage_in=ionos_epotential_mage ) ! transform to pbuf for aurora... @@ -457,6 +479,10 @@ subroutine ionosphere_run1(pbuf2d) call outfld('ltr_efx_phys', pbuf_efx, pcols, lchnk ) call outfld('ltr_kev_phys', pbuf_kev, pcols, lchnk ) end if + if ( ionos_epotential_mage) then + call outfld('mage_efx_phys', pbuf_efx, pcols, lchnk ) + call outfld('mage_kev_phys', pbuf_kev, pcols, lchnk ) + end if end do chnk_loop1 deallocate(prescr_efx, prescr_kev) @@ -1005,6 +1031,9 @@ subroutine ionosphere_init_restart(File) integer :: dimids(3), ndims type(cam_grid_header_info_t) :: info + integer :: mag_grid_id + type(cam_grid_header_info_t) :: mag_grid_info + if (ionos_xport_active) then grid_id = cam_grid_id('physgrid') call cam_grid_write_attr(File, grid_id, info) @@ -1014,11 +1043,21 @@ subroutine ionosphere_init_restart(File) end do ndims = hdimcnt + 1 - call cam_pio_def_dim(File, 'lev', pver, dimids(ndims), & - existOK=.true.) + call cam_pio_def_dim(File, 'lev', pver, dimids(ndims), existOK=.true.) + + ierr = pio_def_var(File, 'Optm1', pio_double, dimids(1:ndims), Optm1_vdesc) + + mag_grid_id = cam_grid_id('gmag_grid') + call cam_grid_write_attr(File, mag_grid_id, mag_grid_info) + + ndims = mag_grid_info%num_hdims() + do i = 1,mag_grid_info%num_hdims() + dimids(i) = mag_grid_info%get_hdimid(i) + end do + + ierr = pio_def_var(File, 'azigm1', pio_double, dimids(1:ndims), azigm1_vdesc) + ierr = pio_def_var(File, 'azigm2', pio_double, dimids(1:ndims), azigm2_vdesc) - ierr = pio_def_var(File, 'Optm1', pio_double, dimids(1:ndims), & - Optm1_vdesc) end if end subroutine ionosphere_init_restart @@ -1030,6 +1069,8 @@ subroutine ionosphere_write_restart(File) use cam_grid_support, only: cam_grid_id, cam_grid_write_var use cam_grid_support, only: cam_grid_get_decomp, cam_grid_dimensions use phys_grid, only: phys_decomp + use edynamo, only: azigm1, azigm2 + use edyn_mpi, only: mlon0,omlon1, mlat0,mlat1 type(file_desc_t), intent(inout) :: File @@ -1039,6 +1080,9 @@ subroutine ionosphere_write_restart(File) integer :: nhdims type(io_desc_t), pointer :: iodesc3d + integer :: mag_grid_id + type(io_desc_t), pointer :: iodesc2d + if (ionos_xport_active) then ! Write grid vars @@ -1051,10 +1095,25 @@ subroutine ionosphere_write_restart(File) dims(1) = pcols dims(2) = pver dims(3) = endchunk - begchunk + 1 - call cam_grid_get_decomp(physgrid, dims(1:3), gdims(1:nhdims), & - pio_double, iodesc3d) + call cam_grid_get_decomp(physgrid, dims(1:3), gdims(1:nhdims), pio_double, iodesc3d) call pio_write_darray(File, Optm1_vdesc, iodesc3d, opmmrtm1_phys, ierr) + + mag_grid_id = cam_grid_id('gmag_grid') + + ! write coords for GoeMag grid + call cam_grid_write_var(File, mag_grid_id) + + ! setup IO decomp + call cam_grid_dimensions(mag_grid_id, gdims(1:2), nhdims) + dims(1) = omlon1 - mlon0 + 1 + dims(2) = mlat1 - mlat0 + 1 + call cam_grid_get_decomp(mag_grid_id, dims(1:2), gdims(1:nhdims), pio_double, iodesc2d) + + ! Write fields on geo-mag grid + call pio_write_darray(File, azigm1_vdesc, iodesc2d, azigm1(mlon0:omlon1,mlat0:mlat1), ierr) + call pio_write_darray(File, azigm2_vdesc, iodesc2d, azigm2(mlon0:omlon1,mlat0:mlat1), ierr) + end if end subroutine ionosphere_write_restart @@ -1064,8 +1123,11 @@ end subroutine ionosphere_write_restart subroutine ionosphere_read_restart(File) use pio, only: io_desc_t, file_desc_t, pio_inq_varid use pio, only: pio_read_darray, pio_double + use pio, only: pio_seterrorhandling, PIO_BCAST_ERROR, PIO_NOERR use cam_grid_support, only: cam_grid_id use cam_grid_support, only: cam_grid_get_decomp, cam_grid_dimensions + use edynamo, only: azigm1, azigm2 + use edyn_mpi, only: mlon0,omlon1, mlat0,mlat1 type(file_desc_t), intent(inout) :: File @@ -1074,6 +1136,9 @@ subroutine ionosphere_read_restart(File) integer :: dims(3), gdims(3) integer :: nhdims type(io_desc_t), pointer :: iodesc3d + type(io_desc_t), pointer :: iodesc2d + integer :: mag_grid_id + integer :: err_handling if (ionos_xport_active) then call ionosphere_alloc() @@ -1085,12 +1150,30 @@ subroutine ionosphere_read_restart(File) dims(1) = pcols dims(2) = pver dims(3) = endchunk - begchunk + 1 - call cam_grid_get_decomp(physgrid, dims(1:3), gdims(1:nhdims), & - pio_double, iodesc3d) + call cam_grid_get_decomp(physgrid, dims(1:3), gdims(1:nhdims), pio_double, iodesc3d) ierr = pio_inq_varid(File, 'Optm1', Optm1_vdesc) call pio_read_darray(File, Optm1_vdesc, iodesc3d, opmmrtm1_phys, ierr) opmmrtm1_initialized = .true. + + mag_grid_id = cam_grid_id('gmag_grid') + call cam_grid_dimensions(mag_grid_id, gdims(1:2), nhdims) + dims(1) = omlon1 - mlon0 + 1 + dims(2) = mlat1 - mlat0 + 1 + call cam_grid_get_decomp(mag_grid_id, dims(1:2), gdims(1:nhdims), pio_double, iodesc2d) + + ! handle errors ourselves + call pio_seterrorhandling(File, PIO_BCAST_ERROR, err_handling) + + ! read vars if available on restart file + ierr = pio_inq_varid(File, 'azigm1', azigm1_vdesc) + if (ierr.eq.PIO_NOERR) call pio_read_darray(File, azigm1_vdesc, iodesc2d, azigm1(mlon0:omlon1,mlat0:mlat1), ierr) + ierr = pio_inq_varid(File, 'azigm2', azigm2_vdesc) + if (ierr.eq.PIO_NOERR) call pio_read_darray(File, azigm2_vdesc, iodesc2d, azigm2(mlon0:omlon1,mlat0:mlat1), ierr) + + ! restore old error handling + call pio_seterrorhandling(File, err_handling) + end if end subroutine ionosphere_read_restart diff --git a/src/ionosphere/waccmx/mage_module.F90 b/src/ionosphere/waccmx/mage_module.F90 new file mode 100644 index 0000000000..fb47d59f92 --- /dev/null +++ b/src/ionosphere/waccmx/mage_module.F90 @@ -0,0 +1,556 @@ +module mage_module + ! + ! Module used to exchange data back and forth with MAGE + ! WACCM-X Receives (POT, mean energy, energy flux) + ! MAGE Receives (Pedersen Conductance, Hall Conductance) + ! + + use shr_kind_mod, only: r8 => shr_kind_r8, cl => shr_kind_cl + use cam_logfile, only: iulog + use spmd_utils, only: masterproc + use edyn_maggrid, only: nmlat, nmlonp1 + use edyn_maggrid, only: gmlon ! magnetic latitudes (nmlat) (radians) + use edyn_maggrid, only: gmlat ! magnetic longtitudes (nmlonp1) (radians) + use edyn_mpi, only: ntask, mytid + use edyn_params, only: pi, dtr, rtd + use edynamo, only: azigm1, azigm2 ! Hall and Ped conductances + use cam_abortutils, only: endrun + + implicit none + + private + public :: mage_init + public :: mage_advance + + integer, parameter :: &! For MPI-based coupling tag + myAppId = 67, &! waccmxID + voltId = 116, &! voltronID + gamId = 45, &! gameraID + rcmId = 34, &! rcmID + mageId = 26, &! mageID + hidraId = 40, &! hidraID + hidraNId = 54, &! hidraNID + hidraSId = 59, &! hidraSID + tiegcmId = 57 ! tiegcmID + + ! Change the following parameters if adding more variables + integer, parameter :: & + nmixinapex = 3, &! APEX Potential, Flux (for auroral bc) + nmixingeo = 0, &! GEO Flux, Energy + nmixoutapex= 2, &! APEX Pedersen, Hall + nmixoutgeo = 0, &! No GEO exports + nhoutvar = 9, &! TN, UN, VN, OMEGA, O2, O1, NO, Z, HE + nhinvar = 0 ! + + integer,dimension(:),allocatable :: IAm + + integer :: & + CplComm, CplCommSize, CplRank, & ! Global coupling communicator + mixCplRank,hidraCplRank, & ! Direct mix coupling communicator + hidraNCplRank,hidraSCplRank ! Direct hidra coupling communicator + + integer :: nlatp2=0, nlonp1=0 ! Needed for geographic conversions later + +contains + + !----------------------------------------------------------------------- + + subroutine mage_advance(iyear, imo, iday, iutsec, & + phihm, mage_efxm, mage_kevm) + + use cam_history_support, only: fillvalue + + ! + ! Read MAGE outputs from mage_ncfile file, returning electric potential, + ! auroral mean energy and energy flux at current date and time, + ! and the data is linearly interpolated to the model time + ! + ! + ! Args: + + integer, intent(in) :: iyear + integer, intent(in) :: imo + integer, intent(in) :: iday + integer, intent(in) :: iutsec + real(r8), intent(out) :: phihm(nmlonp1,nmlat) + real(r8), intent(out) :: mage_efxm(nmlonp1,nmlat) ! on geomag grid + real(r8), intent(out) :: mage_kevm(nmlonp1,nmlat) ! on geomag grid + + if (mytid .eq. 0) then + write(iulog,'(A,I4,1X,I2.2,1X,I2.2,1X,I2,":",I2.2,":",I2.2)') "WCMX Coupling at ", & + iyear, imo, iday, iutsec/3600, mod(iutsec,3600)/60, mod(iutsec,60) + endif + + phihm = 0._r8 + mage_efxm = 0._r8 + mage_kevm = 0._r8 + + call update_mage(phihm, mage_efxm, mage_kevm) + + end subroutine mage_advance + + !----------------------------------------------------------------------- + + subroutine mage_init() + use mpi, only: MPI_DOUBLE_PRECISION, MPI_INTEGER, MPI_ERROR + + integer :: i,ierr + character(len=*), parameter :: preface = 'mage_module::mage_init : WACCMX send to REMIX : ' + + ! Initialize the MPI Coupling interface + + call mp_coupling() + if (mytid .eq. 0) then + i = 0 + if (masterproc) write(iulog,*) "WCMX: Start init_mpi_remix" + + i = i + 1 + call mpi_send(nmlat, 1, MPI_INTEGER, mixCplRank, (myAppId+voltId)*100+i, CplComm, ierr) + if (ierr == MPI_ERROR) then + call endrun(preface//'nmlat') + end if + + i = i + 1 + call mpi_send(nmlonp1, 1, MPI_INTEGER, mixCplRank, (myAppId+voltId)*100+i, CplComm, ierr) + if (ierr == MPI_ERROR) then + call endrun(preface//'nmlonp1') + end if + + i = i + 1 + call mpi_send(gmlat, nmlat, MPI_DOUBLE_PRECISION, mixCplRank, (myAppId+voltId)*100+i, CplComm, ierr) + if (ierr == MPI_ERROR) then + call endrun(preface//'gmlat') + end if + + i = i + 1 + ! Making a note here. For whatever reason, WACCM-X's conductance is flipped 180 degrees compared to TIEGCM. + !I'm going to hack it here so that it's straight up 180 degrees shifted but I do not know why and what the consequences are. + call mpi_send(gmlon, nmlonp1, MPI_DOUBLE_PRECISION, mixCplRank, (myAppId+voltId)*100+i, CplComm, ierr) + if (ierr == MPI_ERROR) then + call endrun(preface//'gmlon, nmlonp1') + end if + + if (masterproc) write(iulog,*) "WCMX: Done init_mpi_remix" + endif + + end subroutine mage_init + + !----------------------------------------------------------------------- + + subroutine mp_coupling() + use mpi, only: mpi_comm_size, mpi_comm_free, MPI_COMM_WORLD, MPI_ERROR, MPI_COMM_NULL, MPI_INTEGER + + integer :: ierr,color,i + integer :: tmpComm + character(len=*), parameter :: preface = 'mage_module::mp_coupling : ' + + ! Create a second communicator to transfer data between TIEGCM and MAGE + ! This communicator only includes the root processes + ! TIEGCM root sends/receives data from/to MAGE root + + color = mageId + tmpComm = MPI_COMM_WORLD + + if (masterproc) write(iulog,*) "WACCMX: Get Coupling Comm ",ntask + + call mp_get_coupling_comm(tmpComm,color, 0, CplComm) + + call mpi_comm_size(CplComm,CplCommSize,ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from mpi_comm_size: ierr=',i4)") ierr + call endrun(preface//'mpi_comm_size') + endif + + ! At most two processes will register in CplComm (TIEGCM root, REMIX root) + ! Therefore CplCommSize is either 1 or 2 + if (CplCommSize == 1) then + + ! Only one process registered in CplComm (TIEGCM root) + ! No coupling will take place, free up the resources + call mpi_comm_free(CplComm,ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from mpi_comm_free: ierr=',i4)") ierr + call endrun(preface//'mpi_comm_free') + endif + CplComm = MPI_COMM_NULL + else + + ! There is another process registered in CplComm (REMIX root) + call mpi_comm_rank(CplComm,CplRank,ierr) + !write(*,*) "W cplrank: ",CplRank + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from mpi_comm_rank: ierr=',i4)") ierr + call endrun(preface//'mpi_comm_rank') + endif + + if (.not.allocated(IAm)) allocate(IAm(CplCommSize)) + + if (mytid .eq. 0) then + IAm(CplRank+1) = myAppId*100+1 + else + IAm(CplRank+1) = myAppId*100 + endif + + do i=1,CplCommSize + call mpi_bcast(IAm(i), 1, MPI_INTEGER, i-1, CplComm, ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from mpi_bcast: ierr=',i4)") ierr + call endrun(preface//'mpi_bcast') + endif + enddo + + do i=1,CplCommSize + ! Assign rank if match + select case (IAm(i)) + case (voltId) + mixCplRank = i-1 + if (mytid .eq. 0) write(iulog,*) "W coupling to remix",mixCplRank + case (gamId) + if (mytid .eq. 0) write(iulog,*) "W not coupling to Gam yet" + case (rcmId) + if (mytid .eq. 0) write(iulog,*) "W not coupling to RCM yet" + case (hidraNId) + hidraNCplRank = i-1 + if (mytid .eq. 0) write(iulog,*) "W coupling to hidraN" + case (hidraSId) + hidraSCplRank = i-1 + if (mytid .eq. 0) write(iulog,*) "W coupling to hidraS" + case (hidraId) + hidraCplRank = i-1 + if (mytid .eq. 0) write(iulog,*) "W coupling to hidra" + case (myAppId,myAppId*100,myAppId*100+1) + if (mytid .eq. 0 .and. IAm(i) .eq. myAppId*100+1) write(iulog,*) "W is W", i-1 + IAm(i) = myAppId + case default + if (mytid .eq. 0) write(iulog,*) "W does not know about", & + " this Coupling ID: ", IAm(i) + end select + enddo + endif + + if (masterproc) write(iulog,'(A,I0,A,I0,A,I0,A)') "W COUPLING to ",CplCommSize, & + " Models on ",CplRank," Rank on ",CplComm," Comm" + + end subroutine mp_coupling + + !----------------------------------------------------------------------- + + subroutine mp_get_coupling_comm(couplingPool,appId,key,coupledComm) + use mpi,only: MPI_comm_rank, MPI_comm_split, MPI_ERROR, MPI_IN_PLACE, MPI_INTEGER, MPI_MAX + + integer, intent(inout) :: couplingPool + integer, intent(in) :: appId, key + integer, intent(inout) :: coupledComm + + integer :: ierr, myRank, appIdCpy + character(len=*), parameter :: preface = 'mage_module::mp_get_coupling_comm : ' + + appIdCpy = appId + + ! mpi_bcast doesn't interact well with intent(in) + + ! tell everyone I'm the broadcasting root + ! broadcast which app I'm creating a communicator with, split with it, and then + ! create a smaller pool that excludes that app + call MPI_comm_rank(couplingPool, myRank, ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from MPI_comm_rank: ierr=',i4)") ierr + call endrun(preface//'MPI_comm_rank') + endif + + call MPI_Allreduce(MPI_IN_PLACE, myRank, 1, MPI_INTEGER, MPI_MAX, couplingPool, ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from MPI_Allreduce: ierr=',i4)") ierr + call endrun(preface//'MPI_Allreduce') + endif + + ! This Bcast is causing a lot of issues. I don't know if this is needed or + ! if it will cause problems for voltron and other models. The behavior here is odd. + call MPI_Bcast(appIdCpy, 1, MPI_INTEGER, myRank, couplingPool, ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from MPI_Bcast: ierr=',i4)") ierr + call endrun(preface//'MPI_Bcast') + endif + + call MPI_comm_split(couplingPool, appIdCpy, key, coupledComm, ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from MPI_comm_split: ierr=',i4)") ierr + call endrun(preface//'MPI_comm_split') + endif + + ! key is never used when making the exclusion pool, 0 is used to preserve order + call MPI_comm_split(couplingPool, myAppId, 0, couplingPool, ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from MPI_comm_split: ierr=',i4)") ierr + call endrun(preface//'MPI_comm_split') + endif + + end subroutine mp_get_coupling_comm + + !----------------------------------------------------------------------- + + subroutine update_mage(phihm, mage_efxm, mage_kevm) + + real(r8), intent(out) :: phihm(nmlonp1,nmlat) + real(r8), intent(out) :: mage_efxm(nmlonp1,nmlat) ! on geomag grid + real(r8), intent(out) :: mage_kevm(nmlonp1,nmlat) ! on geomag grid + + call import_mage(phihm, mage_efxm, mage_kevm) + call export_mage() + + end subroutine update_mage + + !----------------------------------------------------------------------- + + subroutine import_mage(phihm, mage_efxm, mage_kevm) + + real(r8), intent(out) :: phihm(nmlonp1,nmlat) + real(r8), intent(out) :: mage_efxm(nmlonp1,nmlat) ! on geomag grid + real(r8), intent(out) :: mage_kevm(nmlonp1,nmlat) ! on geomag grid + + ! 1. Receive arrays (pot, eng, flx) from REMIX + ! 2. Broadcast to all MPI tasks + ! 3. Set periodic boundaries for gpot, geng, gflx, gpotm + ! 4. Clean up pole values for the dynamo solver + + integer :: i,j + real(r8),dimension(nlatp2,nlonp1,nmixingeo) :: gvar2d + real(r8),dimension(nmlat,nmlonp1,nmixinapex) :: avar2d + + avar2d = 0.0_r8 + gvar2d = 0.0_r8 + + call import_remix(avar2d,gvar2d) + + ! Unit conversion + avar2d(:,:,1) = avar2d(:,:,1)*1e3_r8 ! potential: kV -> V + + ! Process the imported data + + if (nmixoutapex .ne. 0) then + do j=1,nmlat + do i=1,nmlonp1 + phihm(i,j) = avar2d(j,i,1) + !if (avar2d(j,i,2) .ne. avar2d(j,i,2)) write(*,*) "AVAR2?? ",avar2d(j,i,2),j,i + !if (avar2d(j,i,3) .ne. avar2d(j,i,3)) write(*,*) "AVAR3?? ",avar2d(j,i,3),j,i + !if (avar2d(j,i,2) .lt. 0) write(*,*) "AVAR4?? ",avar2d(j,i,2),j,i + !if (avar2d(j,i,3) .lt. 0) write(*,*) "AVAR5?? ",avar2d(j,i,3),j,i + mage_kevm(i,j) = max(avar2d(j,i,2),0.5_r8) ! keV + mage_efxm(i,j) = max(avar2d(j,i,2) * avar2d(j,i,3),0.01_r8)*1.602e-9_r8 ! Convert from keV/cm^s to mW/m^2 + enddo + enddo + endif + + end subroutine import_mage + + !------------------------------------------------------------------- + + subroutine import_remix(avar2d,gvar2d) + use mpi, only: MPI_DOUBLE_PRECISION, MPI_ERROR + + integer :: ierr + real(r8),dimension(:,:,:) :: avar2d,gvar2d + character(len=*), parameter :: preface = 'mage_module::import_remix : ' + + if (nmixingeo .ne. 0) then + call MPI_BCAST(gvar2d, nlatp2*nlonp1*nmixingeo, MPI_DOUBLE_PRECISION, mixCplRank, CplComm, ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from MPI_BCAST: ierr=',i4)") ierr + call endrun(preface//'MPI_BCAST') + endif + endif + if (nmixinapex .ne. 0) then + call MPI_BCAST(avar2d, nmlat*nmlonp1*nmixinapex, MPI_DOUBLE_PRECISION, mixCplRank, CplComm, ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from MPI_BCAST: ierr=',i4)") ierr + call endrun(preface//'MPI_BCAST') + endif + if (masterproc) then + write(iulog,*) "WCMX Done Import2: ",mixCplRank,nmlat*nmlonp1*nmixinapex + write(iulog,*) "WCMX Import Check1: ",minval(avar2d(:,:,1)),maxval(avar2d(:,:,1)) + write(iulog,*) "WCMX Import Check2: ",minval(avar2d(:,:,2)),maxval(avar2d(:,:,2)) + write(iulog,*) "WCMX Import Check3: ",minval(avar2d(:,:,3)),maxval(avar2d(:,:,3)) + endif + endif + + end subroutine import_remix + + !----------------------------------------------------------------------- + + subroutine export_mage() + + logical :: hidra_prep + real(r8),dimension(:,:,:), allocatable :: mixgeoout,mixapexout + integer :: nreq,i + character(len=*), parameter :: preface = 'mage_module::export_mage : ' + + ! Prepare the export data + hidra_prep = .false. + nreq = 0 + !if (mytid == 0) write(iulog,*) "W Starting Export Prep" + + do i=1,CplCommSize + ! Skip Self + if (i == CplRank+1) continue + ! Assign rank if match + select case (IAm(i)) + case (voltId) + call prep_export_remix(mixapexout,mixgeoout) + nreq = nreq + 1 + case (gamId) + !write(*,*) "T not coupling to Gam yet" + case (rcmId) + !write(*,*) "T not coupling to RCM yet" + case (hidraNId) + if (.not.hidra_prep) then + !write(*,*) "Prep H Export" + !call prep_export_hidra(hidraout) + hidra_prep = .true. + endif + nreq = nreq + 1 + case (hidraSId) + if (.not.hidra_prep) then + !write(*,*) "Prep H Export" + !call prep_export_hidra(hidraout) + hidra_prep = .true. + endif + nreq = nreq + 1 + case (hidraId) + if (.not.hidra_prep) then + !call prep_export_hidra(hidraout) + hidra_prep = .true. + endif + nreq = nreq + 1 + case (myAppId) + !write(*,*) "T is T" + case default + if (IAm(i) .eq. 0) cycle + if (mytid == 0) & + write(iulog,*) "W does not know about this Coupling ID: ", IAm(i) + end select + enddo + + ! Allocate request array here + + ! Send the export data + do i=1,CplCommSize + ! Assign rank if match + select case (IAm(i)) + case (voltId) + call export_remix(mixapexout,mixgeoout) + case (gamId) + !write(*,*) "T not coupling to Gam yet" + case (rcmId) + !write(*,*) "T not coupling to RCM yet" + case (hidraNId) + !write(*,*) "Export to HidraN" + !call export_hidra(hidraout,myAppId+hidraNId,hidraNCplRank) + case (hidraSId) + !write(*,*) "Export to HidraS" + !call export_hidra(hidraout,myAppId+hidraSId,hidraSCplRank) + case (hidraId) + !call export_hidra(hidraout,myAppId+hidraId,hidraCplRank) + case (myAppId) + !write(*,*) "T is T" + case default + if (IAm(i) .eq. 0) cycle + if (mytid == 0) & + write(iulog,*) "W does not know about this Coupling ID: ", IAm(i) + end select + enddo + end subroutine export_mage + + !----------------------------------------------------------------------- + + subroutine export_remix(avar2d,gvar2d) + use mpi, only: MPI_DOUBLE_PRECISION, MPI_ERROR + + real(r8),dimension(:,:,:) :: avar2d,gvar2d + integer :: ierr + character(len=*), parameter :: preface = 'mage_module::export_remix : ' + + ! Export the data + if (mytid .eq. 0) then + !write(iulog,*) "WCMX Waiting to export1:",mixCplRank,(myAppId+voltId)*100,nlatp2*nlonp1*nmixoutgeo + if ( nmixoutgeo .ne. 0) then + call mpi_send(gvar2d, nlatp2*nlonp1*nmixoutgeo, MPI_DOUBLE_PRECISION, mixCplRank, & + (myAppId+voltId)*100, CplComm, ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from mpi_send: ierr=',i4)") ierr + call endrun(preface//'mpi_comm_size') + endif + endif + + !write(iulog,*) "WCMX Waiting to export2:",mixCplRank,(myAppId+voltId)*100,nmlat*nmlonp1*nmixoutapex + !write(iulog,*) "WCMX export3: ",shape(avar2d) + if ( nmixoutapex .ne. 0) then + call mpi_send(avar2d, nmlat*nmlonp1*nmixoutapex, MPI_DOUBLE_PRECISION, mixCplRank, & + (myAppId+voltId)*100, CplComm, ierr) + if (ierr == MPI_ERROR) then + write(6,"('>>> Error from mpi_send: ierr=',i4)") ierr + call endrun(preface//'mpi_comm_size') + endif + endif + !write(iulog,*) "WCMX Done to export2:" + + endif + + end subroutine export_remix + !----------------------------------------------------------------------- + !----------------------------------------------------------------------- + subroutine prep_export_remix(avar2d,gvar2d) + use edyn_mpi,only: mlat0,mlat1,mlon0,mlon1, mp_gather_edyn + + ! Local variables .................................................. + integer :: i,j,v + real(r8),dimension(:,:,:),allocatable :: avar2d,gvar2d + + real(r8) :: amsub(mlon0:mlon1,mlat0:mlat1,nmixoutapex) + real(r8) :: amglb(nmlonp1,nmlat,nmixoutapex) + + ! Begin ............................................................ + + ! Prepare data for export + amsub(:,:,1) = azigm2(mlon0:mlon1,mlat0:mlat1) ! Ped + amsub(:,:,2) = azigm1(mlon0:mlon1,mlat0:mlat1) ! Hall + + where(amsub < 0.2_r8) amsub = 0.2_r8 + + call mp_gather_edyn(amsub,mlon0,mlon1,mlat0,mlat1, & + amglb,nmlonp1,nmlat,nmixoutapex) + + + if (mytid .eq. 0) then + !write(iulog,*) "WCMX: Start Preparing Export" + + if (nmixoutapex .ne. 0) then + if (.not. allocated(avar2d)) then + allocate(avar2d(nmlat,nmlonp1,nmixoutapex)) + endif + + do v=1,nmixoutapex + do j = 1,nmlat + do i = 1,nmlonp1 + avar2d(j,i,v) = amglb(i,j,v) + enddo + enddo + enddo + + endif + + if (nmixoutgeo .ne. 0) then + ! Prepare geographic output arrays + if (.not. allocated(gvar2d)) then + allocate(gvar2d(nlatp2,nlonp1,nmixoutgeo)) + endif + + gvar2d(:,:,1) = 0._r8 ! put stuff here + + endif + endif + + !if (mytid .eq. 0) write(iulog,*) "WCMX: Done Preparing Export" + + end subroutine prep_export_remix + +end module mage_module