#define ESMF_CHECK if (ESMF_LogFoundError(rcToCheck=rc, line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! $Id: InjectorMod.F90,v 1.14 2011/06/30 05:58:24 theurich Exp $ ! !------------------------------------------------------------------------- !BOP ! \subsection{Dune Component} ! ! !DESCRIPTION: ! This is a Dune Component. It computes aeolian transport ! ! !EOP module XDuneMod ! ESMF module use ESMF use iso_c_binding implicit none private ! Private data block (This could be s or par) type dunedata sequence real(ESMF_KIND_r8), allocatable :: hfarray(:,:) real(ESMF_KIND_r8), allocatable :: dh_dtfarray(:,:) real(ESMF_KIND_r8), allocatable :: dhfarray(:,:) real(ESMF_KIND_r8), allocatable :: hcarray(:) real(ESMF_KIND_r8), allocatable :: dh_dtcarray(:) real(ESMF_KIND_r8), allocatable :: dhcarray(:) end type dunedata type wrapper sequence type(dunedata), pointer :: ptr end type wrapper ! Define the fortran interface to c non returning functions interface ! Ignore the return value of strncpy -> subroutine ! "restrict" is always assumed if we do not pass a pointer subroutine duneinit() bind(C) import end subroutine duneinit subroutine dunestep() bind(C) import end subroutine dunestep subroutine duneoutput() bind(C) import end subroutine duneoutput subroutine dunesetdtmax(val) bind(C, name="dunesetdtmax") use iso_c_binding ! pff fortran requires use inside each interface definition real(c_double), value :: val end subroutine dunesetdtmax subroutine dunesetdoublepar(string, val) bind(C, name="dunesetdoublepar") use iso_c_binding ! pff fortran requires use inside each interface definition character(kind=c_char), intent(in) :: string(*) ! array of char, remember to end with c_null_char real(c_double), intent(in), value :: val end subroutine dunesetdoublepar subroutine dunegetdouble2dscalar(string, array, nx, ny) bind(C, name="dunegetdouble2dscalar") use iso_c_binding ! pff fortran requires use inside each interface definition character(kind=c_char), intent(in) :: string(*) ! array of char, remember to end with c_null_char real(c_double), dimension(ny*nx), intent(inout) :: array integer(c_int), intent(inout) :: nx, ny end subroutine dunegetdouble2dscalar subroutine dunesetdouble2dscalar(string, array, nx, ny) bind(C, name="dunesetdouble2dscalar") use iso_c_binding ! pff fortran requires use inside each interface definition character(kind=c_char), intent(in) :: string(*) ! array of char, remember to end with c_null_char real(c_double), dimension(ny*nx), intent(inout) :: array integer(c_int), intent(inout) :: nx, ny end subroutine dunesetdouble2dscalar end interface ! Functions interface integer(c_int) function dunegetintpar(string) bind(C, name="dunegetintpar") use iso_c_binding ! pff fortran requires use inside each interface definition character(kind=c_char) :: string(*) ! array of char, remember to end with c_null_char end function dunegetintpar real(c_double) function dunegetdoublepar(string) bind(C, name="dunegetdoublepar") use iso_c_binding ! pff fortran requires use inside each interface definition character(kind=c_char) :: string(*) ! array of char, remember to end with c_null_char end function dunegetdoublepar real(c_double) function dunegetdtmax() bind(C, name="dunegetdtmax") use iso_c_binding ! pff fortran requires use inside each interface definition end function dunegetdtmax real(c_double) function dunegettime() bind(C, name="dunegettime") use iso_c_binding ! pff fortran requires use inside each interface definition end function dunegettime integer(c_int) function dunegetnx() bind(C,name="dunegetnx") use iso_c_binding end function dunegetnx integer(c_int) function dunegetny() bind(C,name="dunegetny") use iso_c_binding end function dunegetny real(c_double) function dunegetdir() bind(C, name="dunegetdir") use iso_c_binding ! pff fortran requires use inside each interface definition end function dunegetdir end interface ! External entry point which will register the Init, Run, and Finalize ! routines for this Component. public XDune_register contains !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ !BOPI ! !IROUTINE: XDune_register - Set the Init, Run, Final routines ! !INTERFACE: subroutine XDune_register(comp, rc) ! ! !ARGUMENTS: type(ESMF_GridComp) :: comp integer, intent(out) :: rc ! ! !DESCRIPTION: ! User-written registration routine. This is the ! only public entry point for this Component. When this is called by ! a higher level component it will register with the Framework the ! subroutines to be called when the Framework needs to Initialize, ! Run, or Finalize this Component. ! ! The arguments are: ! \begin{description} ! \item[comp] ! A Gridded Component. ! \item[rc] ! Return code; equals {\tt ESMF\_SUCCESS} if there are no errors, ! {\tt ESMF\_FAILURE} othewise. ! \end{description} ! !EOPI ! print *, "Registering XDune module" ! Register the callback routines. ! ! This Component has a 1 phase initialization, and a single ! phase run and finalize. call ESMF_GridCompSetEntryPoint(comp, ESMF_METHOD_INITIALIZE, userRoutine=xdune_init, rc=rc) ESMF_CHECK call ESMF_GridCompSetEntryPoint(comp, ESMF_METHOD_RUN, userRoutine=xdune_run, rc=rc) ESMF_CHECK call ESMF_GridCompSetEntryPoint(comp, ESMF_METHOD_FINALIZE, userRoutine=xdune_final, rc=rc) ESMF_CHECK ! print *, "XDune module: Registered Initialize, Run, and Finalize routines" ! print *, "XDune module: Registered Private Data block for Internal State" end subroutine XDune_register !------------------------------------------------------------------------------ !BOPI ! !IROUTINE: User Initialization routine ! !INTERFACE: subroutine xdune_init(gcomp, importState, exportState, clock, rc) ! ! !ARGUMENTS: type(ESMF_GridComp) :: gcomp type(ESMF_Clock) :: clock type(ESMF_State) :: importState, exportState ! Data structures type(ESMF_Field) :: himportfield,hexportfield, dhimportfield, & dhexportfield, zsmeanimportfield type(ESMF_Grid) :: grid type(ESMF_DistGrid) :: distgrid ! For mpi type(ESMF_Array) :: harray, dharray, coordarray integer, intent(out) :: rc integer :: nx, ny, i,j real(c_double) :: xori, yori, dx, dir, pi real(ESMF_KIND_r8), pointer :: coordX(:,:), coordY(:,:), farrayPtr(:,:), fvectorptr(:) ! Computational boundaries (1 number per dimension for coordinates of rectilinear grids) integer :: ubnd(2), lbnd(2) type(ESMF_LocStream) :: locstream real(ESMF_KIND_r8) :: zsmean integer :: count ! no specific kind... !character :: customConv, customPurp !character (len=*), optional :: customAttrList(:) type(wrapper) :: wrap type(dunedata), pointer :: data ! ! !DESCRIPTION: ! User-supplied Initialization routine. Sets up data space, ! and marks which Fields in the export state can be produced ! by this Component. The Coupler will mark which Fields are ! needed by whatever other Component(s) are coupled with this ! Component. The second phase of the Init process will place ! the actual Field data in the export state. ! ! The arguments are: ! \begin{description} ! \item[gcomp] ! A Gridded Component. ! \item[importState] ! State containing the import list. ! \item[exportState] ! State containing the export list. ! \item[clock] ! Clock describing the external time. ! \item[rc] ! Return code; equals {\tt ESMF\_SUCCESS} if there are no errors. ! \end{description} ! !EOPI ! Allocate local data. allocate(data) call duneinit() ! print *, "XDune <----- We should have a dune grid now" ! Now let's define the export state. ! This should expose the data from xdune. rc = ESMF_SUCCESS nx = dunegetnx() ny = dunegetny() xori = dunegetdoublepar('xori' // C_NULL_CHAR) yori = dunegetdoublepar('yori' // C_NULL_CHAR) dx = dunegetdoublepar('dx' // C_NULL_CHAR) dir = dunegetdir() pi = 4.d0*atan(1.d0) ! Define an array, (not sure why we can't use grid here) Initialise to 0.0 allocate(data%dh_dtfarray(nx,ny) ) allocate(data%hfarray(nx,ny) ) allocate(data%dhfarray(nx,ny) ) ! Allocate the c arrays allocate(data%dh_dtcarray(ny*nx) ) allocate(data%hcarray(ny*nx) ) allocate(data%dhcarray(ny*nx) ) ! Default to 0 data%dh_dtfarray = 0 data%hfarray = 0 data%dhfarray = 0 ! C order is different from Fortran, so pass x and y in reversed order call dunegetdouble2dscalar('dh_dt' // C_NULL_CHAR, data%dh_dtcarray, ny, nx) call dunegetdouble2dscalar('h' // C_NULL_CHAR, data%hcarray, ny, nx) ! Convert the shape of the array to fortran data%dh_dtfarray = reshape(data%dh_dtcarray, (/ny, nx/), order=(/2,1/)) data%hfarray = reshape(data%hcarray, (/ny, nx/), order=(/2,1/)) ! Dune grid grid=ESMF_GridCreateNoPeriDim(maxIndex=(/ny,nx/), & coordSys=ESMF_COORDSYS_CART, & ! Cartesian coordinates (x,y not lat,lon) coordDep1=(/1,2/), coordDep2=(/1,2/), & ! Use coordDep1=(/1,2/),coordDep2=(/1,2/) for curvilinear grids. indexflag=ESMF_INDEX_GLOBAL, & ! Needed for xgrid name="Dune grid", rc=rc) ESMF_CHECK ! Allocate coordinates (Assume center locations for staggering of h) call ESMF_GridAddCoord(grid, & staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) ESMF_CHECK ! Get the local grid coordinates for X call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, & staggerloc=ESMF_STAGGERLOC_CENTER, & computationalLBound=lbnd, computationalUBound=ubnd, & farrayPtr=coordX, rc=rc) ESMF_CHECK ! Get the local grid coordinates for Y call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, & staggerloc=ESMF_STAGGERLOC_CENTER, & computationalLBound=lbnd, computationalUBound=ubnd, & farrayPtr=coordY, rc=rc) ESMF_CHECK do i=lbnd(1),ubnd(1) do j=lbnd(2),ubnd(2) coordX(i,j) = (i-1)*dx*cos(dir*2*pi/360)-(j-1)*dx*sin(dir*2*pi/360)+xori coordY(i,j) = (i-1)*dx*sin(dir*2*pi/360)+(j-1)*dx*cos(dir*2*pi/360)+yori end do end do !BOE ! Set the Dune Grid in the Dune Component !BOC call ESMF_GridCompSet(gcomp, grid=grid, rc=rc) ESMF_CHECK !EOC !EOE ! Get the distgrid, required to build the array call ESMF_GridGet(grid=grid, staggerloc=ESMF_STAGGERLOC_CENTER, distgrid=distgrid, rc=rc) ESMF_CHECK ! Create a field bundle with the h array himportfield = ESMF_FieldCreate(grid=grid, & typekind=ESMF_TYPEKIND_r8, & staggerloc=ESMF_STAGGERLOC_CENTER, & indexflag=ESMF_INDEX_GLOBAL, & name="himport", rc=rc) ESMF_CHECK call ESMF_FieldGet(himportfield, localDe=0, farrayPtr=farrayptr, & rc=rc) !totalLBound=ftlb, totalUBound=ftub, totalCount=ftc, farrayptr = 0.0d0 ESMF_CHECK ! Create a field bundle with the h array hexportfield = ESMF_FieldCreate(grid=grid, & typekind=ESMF_TYPEKIND_r8, & staggerloc=ESMF_STAGGERLOC_CENTER, & indexflag=ESMF_INDEX_GLOBAL, & name="hexport", rc=rc) ESMF_CHECK call ESMF_FieldGet(hexportfield, localDe=0, farrayPtr=farrayptr, & rc=rc) !totalLBound=ftlb, totalUBound=ftub, totalCount=ftc, ESMF_CHECK farrayptr = data%hfarray ! Create a field bundle with the h array dhimportfield = ESMF_FieldCreate(grid=grid, & typekind=ESMF_TYPEKIND_r8, & staggerloc=ESMF_STAGGERLOC_CENTER, & indexflag=ESMF_INDEX_GLOBAL, & name="dhimport", rc=rc) ESMF_CHECK call ESMF_FieldGet(dhimportfield, localDe=0, farrayPtr=farrayptr, & rc=rc) !totalLBound=ftlb, totalUBound=ftub, totalCount=ftc, farrayptr = 0.0d0 ESMF_CHECK ! Create a field bundle with the h array dhexportfield = ESMF_FieldCreate(grid=grid, & typekind=ESMF_TYPEKIND_r8, & staggerloc=ESMF_STAGGERLOC_CENTER, & indexflag=ESMF_INDEX_GLOBAL, & name="dhexport", rc=rc) ESMF_CHECK call ESMF_FieldGet(dhexportfield, localDe=0, farrayPtr=farrayptr, & rc=rc) !totalLBound=ftlb, totalUBound=ftub, totalCount=ftc, farrayptr = 0.0d0 ESMF_CHECK ! Creating the LocStream count = 1 locstream = ESMF_LocStreamCreate(name="tide", localCount=count,rc=rc) ! Creating a field for the LocStream zsmean = 0.0d0 zsmeanimportfield = ESMF_FieldCreate(locstream, & ESMF_TYPEKIND_r8, & name="zsmeanimport", rc=rc) ! CHECK ESMF_CHECK ! Temporarily fil it with 0 call ESMF_FieldGet(zsmeanimportfield, localDe=0, farrayPtr=fvectorptr, & rc=rc) farrayptr = 0.0d0 ESMF_CHECK ! Is the state already there? It should have been set by the xdunexbeach component ! Add the field bundle to the import state call ESMF_StateAdd(importstate, (/dhimportfield, himportfield, zsmeanimportfield/), rc=rc) ESMF_CHECK ! Add the field bundle to the epxort state call ESMF_StateAdd(exportstate, (/dhexportfield, hexportfield/), rc=rc) ESMF_CHECK ! Allocate private persistent space (store arrays that are passed to dune in data so they can persist through run and be cleaned in finalize) wrap%ptr => data call ESMF_GridCompSetInternalState(gcomp, wrap, rc) ESMF_CHECK end subroutine xdune_init !------------------------------------------------------------------------------ !BOPI ! !IROUTINE: xdune_run - xdune run routine ! !INTERFACE: subroutine xdune_run(comp, importState, exportState, clock, rc) ! ! !ARGUMENTS: type(ESMF_GridComp) :: comp type(ESMF_State) :: importState, exportState type(ESMF_Clock) :: clock integer, intent(out) :: rc type(ESMF_TimeInterval) :: timeinterval real*8 :: step, tesmf, txbeach, t, told real(c_double) :: xori, yori, dx, dir, pi real(ESMF_KIND_r8), pointer :: coordX(:,:), coordY(:,:), farrayPtr(:,:), fvectorptr(:) type(ESMF_FieldBundle) :: fieldbundle type(ESMF_Field) :: field type(ESMF_Grid) :: grid type(ESMF_DistGrid) :: distgrid ! For mpi type(ESMF_Array) :: array real(ESMF_KIND_r8) :: zsmean ! CHECK DIM real(ESMF_KIND_r8), pointer :: zsmeanimportfdblPtr(:) ! I think this should be a vector... real(ESMF_KIND_r8), allocatable :: dh_dtfarray(:,:), dhfarray(:,:), hfarray(:,:) real(ESMF_KIND_r8), allocatable :: dh_dtcarray(:), dhcarray(:) integer :: ubnd(2), lbnd(2), i type(wrapper) :: wrap type(dunedata), pointer :: data integer :: nx, ny ! ! !DESCRIPTION: ! User-supplied Run routine. Computes the changes to the ! morphology due to wind transport ! ! The arguments are: ! \begin{description} ! \item[comp] ! A Gridded Component. ! \item[importState] ! State containing the import list. ! \item[exportState] ! State containing the export list. ! \item[clock] ! Clock describing the external time. ! \item[rc] ! Return code; equals {\tt ESMF\_SUCCESS} if there are no errors, ! otherwise {\tt ESMF\_FAILURE}. ! \end{description} ! !EOPI rc = ESMF_SUCCESS ! Get the simulation time (seconds since tref) from esmf call ESMF_ClockGet(clock, currSimTime=timeinterval, rc=rc) ! Get the double (assuming it also gives the miliseconds...) ! TODO test with non integer second timestep call ESMF_TimeIntervalGet(timeinterval, s_r8=tesmf, rc=rc) ! Get the timestep from esmf, the coupled timestep call ESMF_ClockGet(clock, timestep=timeinterval, rc=rc) call ESMF_TimeIntervalGet(timeinterval, s_r8=step, rc=rc) ! Get the internal state so we don't have to reallocate variables.... ! Can't use options here call ESMF_GridCompGetInternalState(comp, wrap, rc) data => wrap%ptr ! now we have our variables from init back, we can reuse these.... ! Run until current t from xbeach == current t + timeinterval from esmf ! Set the end time in XBeach (I may have to do an updateoutput times here.... call dunesetdtmax(step) ! Not sure if I have to set a new dtmax if the model only runs until half a timestep. ! TODO: Better build in some assertions before runs nx = dunegetnx() ny = dunegetny() ! Process the import state ! Get a the local array data and set the new values to it. call ESMF_StateGet(importState, itemName="dhimport", field=field, rc=rc) call ESMF_FieldGet(field, farrayPtr=farrayPtr, rc=rc) ! Reshape the dh from XBeach to nx,ny data%dhcarray = reshape(reshape(farrayPtr, (/ny,nx/), order=(/2,1/)), (/ny*nx/)) ! Now put back in h_rotate, so it is rotated back before being stored in h call dunesetdouble2dscalar('dh_rotate' // C_NULL_CHAR, data%dhcarray, nx, ny) ! Retrieve the mean tidal water level call ESMF_StateGet(importState, itemName="zsmeanimport", field=field, rc=rc) call ESMF_FieldGet(field, farrayPtr=zsmeanimportfdblPtr, rc=rc) ! Set the tidal water level in Dune ! I expect that there is only 1 value here.... call dunesetdoublepar('tide' // C_NULL_CHAR, zsmeanimportfdblPtr(1)) allocate(dh_dtfarray(nx,ny) ) ! Allocate in c order allocate(dh_dtcarray(ny*nx) ) allocate(dhfarray(nx,ny) ) dhfarray = 0 t = dunegettime() do while (t < tesmf+step) told = dunegettime() call dunestep() t = dunegettime() call dunegetdouble2dscalar('dh_dt_rotate' // C_NULL_CHAR, dh_dtcarray, nx, ny) dh_dtfarray = reshape(dh_dtcarray, (/ nx, ny /), order = (/ 2,1/)) write(*,*) minval(dh_dtfarray), maxval(dh_dtfarray) dhfarray = dhfarray + (dh_dtfarray * (t - told)) end do call duneoutput() ! Get the dh array from the state ! Get a the local array data and set the new values to it. call ESMF_StateGet(exportState, itemName="dhexport", field=field, rc=rc) call ESMF_FieldGet(field, farrayPtr=farrayPtr, rc=rc) data%dhfarray = dhfarray ! Make sure the pointer persists past this subroutine call farrayPtr = data%dhfarray ! Get a the local array data and set the new values to it. call ESMF_StateGet(exportState, itemName="hexport", field=field, rc=rc) call ESMF_FieldGet(field, farrayPtr=farrayPtr, rc=rc) call dunegetdouble2dscalar('h' // C_NULL_CHAR, data%hcarray, nx, ny) data%hfarray = reshape(data%hcarray, (/nx,ny/), order=(/2,1/)) farrayPtr = data%hfarray dir = dunegetdir() pi = 4.d0*atan(1.d0) !GridCompGet call ESMF_GridCompGet(comp, grid=grid, rc=rc) ESMF_CHECK ! call ESMF_GridGet(grid=grid, rc=rc) ! ESMF_CHECK ! Get the local grid coordinates for X call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, & staggerloc=ESMF_STAGGERLOC_CENTER, & computationalLBound=lbnd, computationalUBound=ubnd, & farrayPtr=coordX, rc=rc) ESMF_CHECK do i=lbnd(1),ubnd(1) coordX(i,:) = (i-1)*dx*cos(dir*2*pi/360)+xori end do ! Get the local grid coordinates for Y call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, & staggerloc=ESMF_STAGGERLOC_CENTER, & computationalLBound=lbnd, computationalUBound=ubnd, & farrayPtr=coordY, rc=rc) ESMF_CHECK ! Set the coordinates for y do i=lbnd(2),ubnd(2) coordY(:,i) = (i-1)*dx*sin(dir*2*pi/360)+yori end do call ESMF_GridCompSetInternalState(comp, wrap, rc) write(*,*) 'XDune ends at ', t deallocate(dhfarray) deallocate(dh_dtfarray) deallocate(dh_dtcarray) end subroutine xdune_run !------------------------------------------------------------------------------ !BOPI ! !IROUTINE: xdune_final - finalize routine ! !INTERFACE: subroutine xdune_final(comp, importState, exportState, clock, rc) ! ! !ARGUMENTS: type(ESMF_GridComp) :: comp type(ESMF_State) :: importState, exportState type(ESMF_Clock) :: clock integer, intent(out) :: rc ! ! !DESCRIPTION: ! User-supplied finalize routine. Release space allocated ! by this component. ! ! The arguments are: ! \begin{description} ! \item[comp] ! A Gridded Component. ! \item[importState] ! State containing the import list. ! \item[exportState] ! State containing the export list. ! \item[clock] ! Clock describing the external time. ! \item[rc] ! Return code; equals {\tt ESMF\_SUCCESS} if there are no errors, ! otherwise {\tt ESMF\_FAILURE}. ! \end{description} ! !EOPI rc = ESMF_SUCCESS call ESMF_StateDestroy(exportState, rc=rc) end subroutine xdune_final end module XDuneMod