diff --git a/AdvCore_GridCompMod.F90 b/AdvCore_GridCompMod.F90 index ea45be0..c4ce159 100755 --- a/AdvCore_GridCompMod.F90 +++ b/AdvCore_GridCompMod.F90 @@ -56,7 +56,7 @@ module AdvCore_GridCompMod ! !USES: use ESMF - use MAPL_Mod + use MAPL use m_set_eta, only: set_eta use fv_arrays_mod, only: fv_atmos_type, FVPRC, REAL4, REAL8 use fms_mod, only: fms_init, set_domain, nullify_domain @@ -67,6 +67,7 @@ module AdvCore_GridCompMod USE FV_StateMod, only: AdvCoreTracers => T_TRACERS USE FV_StateMod, only: FV_Atm + use CubeGridPrototype, only: register_grid_and_regridders implicit none private @@ -274,6 +275,7 @@ subroutine SetServices(GC, rc) if (.NOT. FV3_DynCoreIsRunning) then ! Make sure FV3 is setup ! ----------------------- + call register_grid_and_regridders() call fv_init1(FV_Atm, dt, grids_on_my_pe, p_split) ! Get Domain decomposition !------------------------- @@ -316,7 +318,6 @@ subroutine SetServices(GC, rc) !-------------------------------------------------- if (.NOT. FV3_DynCoreIsRunning) then call fv_init2(FV_Atm, dt, grids_on_my_pe, p_split) - call register_grid_and_regridders() end if ! Ending with a Generic SetServices call is a MAPL requirement @@ -935,31 +936,4 @@ subroutine global_integral (QG,Q,PLE,IM,JM,KM,NQ) end subroutine global_integral -subroutine register_grid_and_regridders() - use MAPL_GridManagerMod, only: grid_manager - use CubedSphereGridFactoryMod, only: CubedSphereGridFactory - use MAPL_RegridderManagerMod, only: regridder_manager - use MAPL_RegridderSpecMod, only: REGRID_METHOD_BILINEAR - use LatLonToCubeRegridderMod - use CubeToLatLonRegridderMod - use CubeToCubeRegridderMod - - type (CubedSphereGridFactory) :: factory - - type (CubeToLatLonRegridder) :: cube_to_latlon_prototype - type (LatLonToCubeRegridder) :: latlon_to_cube_prototype - type (CubeToCubeRegridder) :: cube_to_cube_prototype - - call grid_manager%add_prototype('Cubed-Sphere',factory) - associate (method => REGRID_METHOD_BILINEAR, mgr => regridder_manager) - call mgr%add_prototype('Cubed-Sphere', 'LatLon', method, cube_to_latlon_prototype) - call mgr%add_prototype('LatLon', 'Cubed-Sphere', method, latlon_to_cube_prototype) - call mgr%add_prototype('Cubed-Sphere', 'Cubed-Sphere', method, cube_to_cube_prototype) - end associate - - end subroutine register_grid_and_regridders - -!EOC -!------------------------------------------------------------------------------ - end module AdvCore_GridCompMod diff --git a/AppGridCreate.F90 b/AppGridCreate.F90 index 69f33ae..2d9551a 100644 --- a/AppGridCreate.F90 +++ b/AppGridCreate.F90 @@ -3,302 +3,299 @@ subroutine AppCSEdgeCreateF(IM_WORLD, LonEdge,LatEdge, LonCenter, LatCenter, rc) #include "MAPL_Generic.h" - use ESMF - use MAPL_BaseMod - use MAPL_GenericMod - use MAPL_ConstantsMod, only : pi=> MAPL_PI_R8 - use fv_arrays_mod, only: REAL4, REAL8, R_GRID - use fv_grid_utils_mod, only: gnomonic_grids, cell_center2, direct_transform - use fv_grid_tools_mod, only: mirror_grid - use FV_StateMod, only: FV_Atm - implicit none - -! !ARGUMENTS: - integer, intent(IN) :: IM_WORLD - integer, optional, intent(OUT) :: rc - real(ESMF_KIND_R8), intent(inout) :: LonEdge(IM_World+1,IM_World+1,6) - real(ESMF_KIND_R8), intent(inout) :: LatEdge(IM_World+1,IM_World+1,6) - real(ESMF_KIND_R8), optional, intent(inout) :: LonCenter(IM_World,IM_World) - real(ESMF_KIND_R8), optional, intent(inout) :: LatCenter(IM_World,IM_World) - -! ErrLog variables -!----------------- - - integer :: STATUS - character(len=ESMF_MAXSTR), parameter :: Iam="AppCSEdgeCreateF" - -! Local variables -!----------------- + use ESMF + use MAPL + use MAPL_ConstantsMod, only : pi=> MAPL_PI_R8 + use fv_arrays_mod, only: REAL4, REAL8, R_GRID + use fv_grid_utils_mod, only: gnomonic_grids, cell_center2, direct_transform + use fv_grid_tools_mod, only: mirror_grid + use FV_StateMod, only: FV_Atm + implicit none + + ! !ARGUMENTS: + integer, intent(IN) :: IM_WORLD + integer, optional, intent(OUT) :: rc + real(ESMF_KIND_R8), intent(inout) :: LonEdge(IM_World+1,IM_World+1,6) + real(ESMF_KIND_R8), intent(inout) :: LatEdge(IM_World+1,IM_World+1,6) + real(ESMF_KIND_R8), optional, intent(inout) :: LonCenter(IM_World,IM_World) + real(ESMF_KIND_R8), optional, intent(inout) :: LatCenter(IM_World,IM_World) + + ! ErrLog variables + !----------------- + + integer :: STATUS + character(len=ESMF_MAXSTR), parameter :: Iam="AppCSEdgeCreateF" + + ! Local variables + !----------------- #ifdef EIGHT_BYTE - integer, parameter:: f_p = selected_real_kind(15) ! same as 12 on Altix + integer, parameter:: f_p = selected_real_kind(15) ! same as 12 on Altix #else -! Higher precisions for grid geometrical factors: - integer, parameter:: f_p = selected_real_kind(20) + ! Higher precisions for grid geometrical factors: + integer, parameter:: f_p = selected_real_kind(20) #endif - integer, parameter :: grid_type = 0 - integer :: npts - integer :: ntiles=6 - integer :: ndims=2 - integer :: I, J, N - integer :: IG, JG - real(ESMF_KIND_R8) :: alocs(2) - - real(R_GRID), allocatable :: grid_global(:,:,:,:) - - integer :: L - - npts = IM_World - allocate( grid_global(npts+1,npts+1,ndims,ntiles) ) - call gnomonic_grids(grid_type, npts, grid_global(:,:,1,1), grid_global(:,:,2,1)) -! mirror_grid assumes that the tile=1 is centered on equator and greenwich meridian Lon[-pi,pi] - call mirror_grid(grid_global, 0, npts+1, npts+1, 2, 6) - - if (allocated(FV_Atm)) then - if (.not. FV_Atm(1)%flagstruct%do_schmidt) & - grid_global(:,:,1,:) = grid_global(:,:,1,:) - PI/FV_Atm(1)%flagstruct%shift_fac - else - grid_global(:,:,1,:) = grid_global(:,:,1,:) - PI/18.0 - end if - - where(grid_global(:,:,1,:) < 0.) & - grid_global(:,:,1,:) = grid_global(:,:,1,:) + 2.* PI - - ! Keep Equator and Greenwich exact - !--------------------------------- - - where(abs(grid_global(:,:,:,1)) < 1.d-10) grid_global(:,:,:,1) = 0.0 - -!--------------------------------- -! Clean Up Corners -!--------------------------------- - grid_global( 1,1:npts+1,:,2)=grid_global(npts+1,1:npts+1,:,1) - grid_global( 1,1:npts+1,:,3)=grid_global(npts+1:1:-1,npts+1,:,1) - grid_global(1:npts+1,npts+1,:,5)=grid_global(1,npts+1:1:-1,:,1) - grid_global(1:npts+1,npts+1,:,6)=grid_global(1:npts+1,1,:,1) - grid_global(1:npts+1, 1,:,3)=grid_global(1:npts+1,npts+1,:,2) - grid_global(1:npts+1, 1,:,4)=grid_global(npts+1,npts+1:1:-1,:,2) - grid_global(npts+1,1:npts+1,:,6)=grid_global(npts+1:1:-1,1,:,2) - grid_global( 1,1:npts+1,:,4)=grid_global(npts+1,1:npts+1,:,3) - grid_global( 1,1:npts+1,:,5)=grid_global(npts+1:1:-1,npts+1,:,3) - grid_global(npts+1,1:npts+1,:,3)=grid_global(1,1:npts+1,:,4) - grid_global(1:npts+1, 1,:,5)=grid_global(1:npts+1,npts+1,:,4) - grid_global(1:npts+1, 1,:,6)=grid_global(npts+1,npts+1:1:-1,:,4) - grid_global( 1,1:npts+1,:,6)=grid_global(npts+1,1:npts+1,:,5) + integer, parameter :: grid_type = 0 + integer :: npts + integer :: ntiles=6 + integer :: ndims=2 + integer :: I, J, N + integer :: IG, JG + real(ESMF_KIND_R8) :: alocs(2) + + real(R_GRID), allocatable :: grid_global(:,:,:,:) + + integer :: L + + npts = IM_World + allocate( grid_global(npts+1,npts+1,ndims,ntiles) ) + call gnomonic_grids(grid_type, npts, grid_global(:,:,1,1), grid_global(:,:,2,1)) + ! mirror_grid assumes that the tile=1 is centered on equator and greenwich meridian Lon[-pi,pi] + call mirror_grid(grid_global, 0, npts+1, npts+1, 2, 6) + + if (allocated(FV_Atm)) then + if (.not. FV_Atm(1)%flagstruct%do_schmidt) & + grid_global(:,:,1,:) = grid_global(:,:,1,:) - PI/FV_Atm(1)%flagstruct%shift_fac + else + grid_global(:,:,1,:) = grid_global(:,:,1,:) - PI/18.0 + end if + + where(grid_global(:,:,1,:) < 0.) & + grid_global(:,:,1,:) = grid_global(:,:,1,:) + 2.* PI + + ! Keep Equator and Greenwich exact + !--------------------------------- + + where(abs(grid_global(:,:,:,1)) < 1.d-10) grid_global(:,:,:,1) = 0.0 + + !--------------------------------- + ! Clean Up Corners + !--------------------------------- + grid_global( 1,1:npts+1,:,2)=grid_global(npts+1,1:npts+1,:,1) + grid_global( 1,1:npts+1,:,3)=grid_global(npts+1:1:-1,npts+1,:,1) + grid_global(1:npts+1,npts+1,:,5)=grid_global(1,npts+1:1:-1,:,1) + grid_global(1:npts+1,npts+1,:,6)=grid_global(1:npts+1,1,:,1) + grid_global(1:npts+1, 1,:,3)=grid_global(1:npts+1,npts+1,:,2) + grid_global(1:npts+1, 1,:,4)=grid_global(npts+1,npts+1:1:-1,:,2) + grid_global(npts+1,1:npts+1,:,6)=grid_global(npts+1:1:-1,1,:,2) + grid_global( 1,1:npts+1,:,4)=grid_global(npts+1,1:npts+1,:,3) + grid_global( 1,1:npts+1,:,5)=grid_global(npts+1:1:-1,npts+1,:,3) + grid_global(npts+1,1:npts+1,:,3)=grid_global(1,1:npts+1,:,4) + grid_global(1:npts+1, 1,:,5)=grid_global(1:npts+1,npts+1,:,4) + grid_global(1:npts+1, 1,:,6)=grid_global(npts+1,npts+1:1:-1,:,4) + grid_global( 1,1:npts+1,:,6)=grid_global(npts+1,1:npts+1,:,5) !------------------------ ! Schmidt transformation: !------------------------ if (allocated(FV_Atm)) then - if ( FV_Atm(1)%flagstruct%do_schmidt ) then - do n=1,ntiles - call direct_transform(FV_Atm(1)%flagstruct%stretch_fac, 1, npts+1, 1, npts+1, FV_Atm(1)%flagstruct%target_lon, FV_Atm(1)%flagstruct%target_lat, & - n, grid_global(1:npts+1,1:npts+1,1,n), grid_global(1:npts+1,1:npts+1,2,n)) - enddo - endif + if ( FV_Atm(1)%flagstruct%do_schmidt ) then + do n=1,ntiles + call direct_transform(FV_Atm(1)%flagstruct%stretch_fac, 1, npts+1, 1, npts+1, FV_Atm(1)%flagstruct%target_lon, FV_Atm(1)%flagstruct%target_lat, & + n, grid_global(1:npts+1,1:npts+1,1,n), grid_global(1:npts+1,1:npts+1,2,n)) + enddo + endif end if - if (present(LonCenter) .and. present(LatCenter)) then - do n=1,ntiles - do j=1,npts - do i=1,npts - call cell_center2(grid_global(i,j, 1:2,n), grid_global(i+1,j, 1:2,n), & - grid_global(i,j+1,1:2,n), grid_global(i+1,j+1,1:2,n), & - alocs) - jg = (n-1)*npts + j - LonCenter(i,jg) = alocs(1) - LatCenter(i,jg) = alocs(2) - enddo - enddo - enddo - end if + if (present(LonCenter) .and. present(LatCenter)) then + do n=1,ntiles + do j=1,npts + do i=1,npts + call cell_center2(grid_global(i,j, 1:2,n), grid_global(i+1,j, 1:2,n), & + grid_global(i,j+1,1:2,n), grid_global(i+1,j+1,1:2,n), & + alocs) + jg = (n-1)*npts + j + LonCenter(i,jg) = alocs(1) + LatCenter(i,jg) = alocs(2) + enddo + enddo + enddo + end if - LonEdge = grid_global(:,:,1,:) - LatEdge = grid_global(:,:,2,:) + LonEdge = grid_global(:,:,1,:) + LatEdge = grid_global(:,:,2,:) - deallocate( grid_global ) + deallocate( grid_global ) - RETURN_(ESMF_SUCCESS) - end subroutine AppCSEdgeCreateF + RETURN_(ESMF_SUCCESS) +end subroutine AppCSEdgeCreateF !!!!!!!!!!!!!!!%%%%%%%%%%%%%%%%%%%%%%%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function AppGridCreateF(IM_WORLD, JM_WORLD, LM, NX, NY, rc) result(esmfgrid) #include "MAPL_Generic.h" -#define DEALLOCGLOB_(A) if(associated(A)) then;A=0;call MAPL_DeAllocNodeArray(A,rc=STATUS);if(STATUS==MAPL_NoShm) deallocate(A,stat=STATUS);NULLIFY(A);endif - - use ESMF - use MAPL_Mod - use MAPL_BaseMod - use MAPL_GenericMod - use MAPL_ConstantsMod, only : pi=> MAPL_PI_R8 - use MAPL_ShmemMod - - use fv_arrays_mod, only: REAL4, REAL8, R_GRID - use fv_grid_utils_mod, only: gnomonic_grids, cell_center2, direct_transform - use fv_grid_tools_mod, only: mirror_grid - use FV_StateMod, only: FV_Atm - - implicit none - -! !ARGUMENTS: - integer, intent(IN) :: IM_WORLD, JM_WORLD, LM - integer, intent(IN) :: NX, NY - integer, optional, intent(OUT) :: rc - type (ESMF_Grid) :: esmfgrid - -! ErrLog variables -!----------------- - - integer :: STATUS - character(len=ESMF_MAXSTR), parameter :: Iam="AppGridCreateF" - -! Local variables -!----------------- - - type(ESMF_DistGrid) :: distGrid - integer :: DECOUNT(2) - integer :: GLOBAL_DE_START(2) - integer :: NX2, NY2 - integer :: NPES, NPES_X, NPES_Y - integer :: NPX, NPY - integer :: isg, ieg - integer :: jsg, jeg - integer :: is, ie - integer :: js, je - integer :: myTile - integer :: npts - integer :: ntiles, grid_type - integer :: ndims=2 - integer :: I, J, N - integer :: IG, JG - - real(REAL8) :: deglon=0.0 - type (ESMF_Array), pointer :: coords(:) - real(REAL8), pointer :: lons(:,:) - real(REAL8), pointer :: lats(:,:) - type (ESMF_Array), target :: tarray(2) - real(REAL8) :: alocs(2), dx, dy - - integer :: L - integer, allocatable :: IMS(:), JMS(:) - character(len=ESMF_MAXSTR) :: gridname, FMT, FMTIM, FMTJM - type (ESMF_VM) :: VM - type (ESMF_Config) :: cf - real(REAL8), allocatable :: gridCornerLats(:) - real(REAL8), allocatable :: gridCornerLons(:) - integer :: idx - real(R_GRID), pointer :: grid_global(:,:,:,:) => null() - logical :: AppGridStandAlone - -! We need the VM to create the grid ?? -!------------------------------------ - - AppGridStandAlone = .false. - call ESMF_VMGetCurrent(vm, rc=STATUS) - VERIFY_(STATUS) - -! Check FV3 grid_type -! ------------------- - cf = ESMF_ConfigCreate(rc=rc) - call ESMF_ConfigLoadFile( cf, 'fvcore_layout.rc', rc = rc ) - call ESMF_ConfigGetAttribute( cf, ntiles, label = 'ntiles:', default=6, rc = rc ) - call ESMF_ConfigGetAttribute( cf, grid_type, label = 'grid_type:', default=0, rc = rc ) - call ESMF_ConfigDestroy(cf, rc=rc) - -! Give the IMS and JMS the MAPL default distribution -! -------------------------------------------------- - - allocate( IMS(0:NX-1) ) - allocate( JMS(0:NY-1) ) - - call GET_INT_FORMAT(IM_WORLD, FMTIM) - call GET_INT_FORMAT(JM_WORLD, FMTJM) - FMT = '(A,' // trim(FMTIM) //',A,' // trim(FMTJM) // ',A)' - if (grid_type <= 3) then - ! Cubed-Sphere - call DecomposeDim ( IM_WORLD , IMS , NX ) - call DecomposeDim ( JM_WORLD/6, JMS(0:NY/6 -1) , NY/6 ) - do n=2,6 - JMS((n-1)*NY/6 : n*NY/6 -1) = JMS(0:NY/6 -1) - enddo - write(gridname,trim(FMT)) 'PE',IM_WORLD,'x',JM_WORLD,'-CF' - else - ! Doubly Periodic Cartesian - call DecomposeDim ( IM_WORLD , IMS , NX ) - call DecomposeDim ( JM_WORLD , JMS , NY ) - write(gridname,trim(FMT)) 'PE',IM_WORLD,'x',JM_WORLD,'-DP' - endif - -! We should have a much simpler create with the next ESMF grid design -!-------------------------------------------------------------------- - - esmfgrid = ESMF_GridCreate( & - name=gridname, & - countsPerDEDim1=IMS, & - countsPerDEDim2=JMS, & - indexFlag = ESMF_INDEX_DELOCAL,& - gridEdgeLWidth = (/0,0/), & - gridEdgeUWidth = (/0,0/), & - coordDep1 = (/1,2/), & - coordDep2 = (/1,2/), & - rc=status) - VERIFY_(STATUS) - -! Allocate coords at default stagger location - call ESMF_GridAddCoord(esmfgrid, rc=status) - VERIFY_(STATUS) - - call ESMF_AttributeSet(esmfgrid, name='GRID_LM', value=LM, rc=status) - VERIFY_(STATUS) - - if (grid_type <= 3) then - call ESMF_AttributeSet(esmfgrid, 'GridType', 'Cubed-Sphere', rc=STATUS) - VERIFY_(STATUS) - else - call ESMF_AttributeSet(esmfgrid, 'GridType', 'Doubly-Periodic', rc=STATUS) - VERIFY_(STATUS) - endif - -! ------------- - - NPES_X = size(IMS) - NPES_Y = size(JMS) - NPES = NPES_X+NPES_Y +#define DEALLOCGLOB_(A) if(associated(A))then;A=0;if(MAPL_ShmInitialized)then; call MAPL_DeAllocNodeArray(A,rc=STATUS);else; deallocate(A,stat=STATUS);endif;VERIFY_(STATUS);NULLIFY(A);endif + + use ESMF + use MAPL, pi=> MAPL_PI_R8 + + use fv_arrays_mod, only: REAL4, REAL8, R_GRID + use fv_grid_utils_mod, only: gnomonic_grids, cell_center2, direct_transform + use fv_grid_tools_mod, only: mirror_grid + use FV_StateMod, only: FV_Atm + + implicit none + + ! !ARGUMENTS: + integer, intent(IN) :: IM_WORLD, JM_WORLD, LM + integer, intent(IN) :: NX, NY + integer, optional, intent(OUT) :: rc + type (ESMF_Grid) :: esmfgrid + + ! ErrLog variables + !----------------- + + integer :: STATUS + character(len=ESMF_MAXSTR), parameter :: Iam="AppGridCreateF" + + ! Local variables + !----------------- + + type(ESMF_DistGrid) :: distGrid + integer :: DECOUNT(2) + integer :: GLOBAL_DE_START(2) + integer :: NX2, NY2 + integer :: NPES, NPES_X, NPES_Y + integer :: NPX, NPY + integer :: isg, ieg + integer :: jsg, jeg + integer :: is, ie + integer :: js, je + integer :: myTile + integer :: npts + integer :: ntiles, grid_type + integer :: ndims=2 + integer :: I, J, N + integer :: IG, JG + + real(REAL8) :: deglon=0.0 + type (ESMF_Array), pointer :: coords(:) + real(REAL8), pointer :: lons(:,:) + real(REAL8), pointer :: lats(:,:) + type (ESMF_Array), target :: tarray(2) + real(REAL8) :: alocs(2), dx, dy + + integer :: L + integer, allocatable :: IMS(:), JMS(:) + character(len=ESMF_MAXSTR) :: gridname, FMT, FMTIM, FMTJM + type (ESMF_VM) :: VM + type (ESMF_Config) :: cf + real(REAL8), allocatable :: gridCornerLats(:) + real(REAL8), allocatable :: gridCornerLons(:) + integer :: idx + real(R_GRID), pointer :: grid_global(:,:,:,:) => null() + logical :: AppGridStandAlone + + ! We need the VM to create the grid ?? + !------------------------------------ + + AppGridStandAlone = .false. + call ESMF_VMGetCurrent(vm, rc=STATUS) + VERIFY_(STATUS) + + ! Check FV3 grid_type + ! ------------------- + cf = ESMF_ConfigCreate(rc=rc) + call ESMF_ConfigLoadFile( cf, 'fvcore_layout.rc', rc = rc ) + call ESMF_ConfigGetAttribute( cf, ntiles, label = 'ntiles:', default=6, rc = rc ) + call ESMF_ConfigGetAttribute( cf, grid_type, label = 'grid_type:', default=0, rc = rc ) + call ESMF_ConfigDestroy(cf, rc=rc) + + ! Give the IMS and JMS the MAPL default distribution + ! -------------------------------------------------- + + allocate( IMS(0:NX-1) ) + allocate( JMS(0:NY-1) ) + + call GET_INT_FORMAT(IM_WORLD, FMTIM) + call GET_INT_FORMAT(JM_WORLD, FMTJM) + FMT = '(A,' // trim(FMTIM) //',A,' // trim(FMTJM) // ',A)' + if (grid_type <= 3) then + ! Cubed-Sphere + call DecomposeDim ( IM_WORLD , IMS , NX ) + call DecomposeDim ( JM_WORLD/6, JMS(0:NY/6 -1) , NY/6 ) + do n=2,6 + JMS((n-1)*NY/6 : n*NY/6 -1) = JMS(0:NY/6 -1) + enddo + write(gridname,trim(FMT)) 'PE',IM_WORLD,'x',JM_WORLD,'-CF' + else + ! Doubly Periodic Cartesian + call DecomposeDim ( IM_WORLD , IMS , NX ) + call DecomposeDim ( JM_WORLD , JMS , NY ) + write(gridname,trim(FMT)) 'PE',IM_WORLD,'x',JM_WORLD,'-DP' + endif + + ! We should have a much simpler create with the next ESMF grid design + !-------------------------------------------------------------------- + + esmfgrid = ESMF_GridCreate( & + name=gridname, & + countsPerDEDim1=IMS, & + countsPerDEDim2=JMS, & + indexFlag = ESMF_INDEX_DELOCAL,& + gridEdgeLWidth = (/0,0/), & + gridEdgeUWidth = (/0,0/), & + coordDep1 = (/1,2/), & + coordDep2 = (/1,2/), & + rc=status) + VERIFY_(STATUS) + + ! Allocate coords at default stagger location + call ESMF_GridAddCoord(esmfgrid, rc=status) + VERIFY_(STATUS) + + call ESMF_AttributeSet(esmfgrid, name='GRID_LM', value=LM, rc=status) + VERIFY_(STATUS) + + if (grid_type <= 3) then + call ESMF_AttributeSet(esmfgrid, 'GridType', 'Cubed-Sphere', rc=STATUS) + VERIFY_(STATUS) + else + call ESMF_AttributeSet(esmfgrid, 'GridType', 'Doubly-Periodic', rc=STATUS) + VERIFY_(STATUS) + endif + + ! ------------- + + NPES_X = size(IMS) + NPES_Y = size(JMS) + NPES = NPES_X+NPES_Y call MAPL_GRID_INTERIOR(esmfgrid,isg,ieg,jsg,jeg) - npx = IM_WORLD - npy = JM_WORLD - myTile = jsg/(npy/ntiles) - - is = isg - ie = ieg - js = jsg - myTile*(npy/ntiles) - je = jeg - myTile*(npy/ntiles) - - npts = (npy/ntiles) - if (npts /= npx) then - print*, 'Error npts /= npx', npts, npx - STATUS=1 - endif - VERIFY_(STATUS) - - if (.not.allocated(FV_Atm)) then - - AppGridStandAlone = .true. - call MAPL_AllocNodeArray(grid_global,Shp=(/npts+1,npts+1,ndims,ntiles/),rc=status) - if (status == MAPL_NoShm) then - allocate( grid_global(npts+1,npts+1,ndims,ntiles) ) - end if - - if (grid_type <= 3) then - ! Cubed-Sphere - call gnomonic_grids(grid_type, npts, grid_global(:,:,1,1), grid_global(:,:,2,1)) - ! mirror_grid assumes that the tile=1 is centered on equator and greenwich meridian Lon[-pi,pi] - call mirror_grid(grid_global, 0, npts+1, npts+1, 2, 6) + npx = IM_WORLD + npy = JM_WORLD + myTile = jsg/(npy/ntiles) + + is = isg + ie = ieg + js = jsg - myTile*(npy/ntiles) + je = jeg - myTile*(npy/ntiles) + + npts = (npy/ntiles) + if (npts /= npx) then + print*, 'Error npts /= npx', npts, npx + STATUS=1 + endif + VERIFY_(STATUS) + + if (.not.allocated(FV_Atm)) then + + AppGridStandAlone = .true. + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(grid_global,Shp=(/npts+1,npts+1,ndims,ntiles/),rc=status) + VERIFY_(STATUS) + else + allocate( grid_global(npts+1,npts+1,ndims,ntiles) ) + end if + + if (grid_type <= 3) then + ! Cubed-Sphere + call gnomonic_grids(grid_type, npts, grid_global(:,:,1,1), grid_global(:,:,2,1)) + ! mirror_grid assumes that the tile=1 is centered on equator and greenwich meridian Lon[-pi,pi] + call mirror_grid(grid_global, 0, npts+1, npts+1, 2, 6) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! ! MAT BMA You cannot refer to FV_Atm here because line 286 says this ! + ! ! MAT BMA You cannot refer to FV_Atm here because line 288 says this ! ! ! is only for FV_Atm not allocated ! ! ! ! if (.not. FV_Atm(1)%flagstruct%do_schmidt) & ! @@ -307,33 +304,33 @@ function AppGridCreateF(IM_WORLD, JM_WORLD, LM, NX, NY, rc) result(esmfgrid) grid_global(:,:,1,:) = grid_global(:,:,1,:) - PI/18.0 - where(grid_global(:,:,1,:) < 0.) & - grid_global(:,:,1,:) = grid_global(:,:,1,:) + 2.* PI - - ! Keep Equator and Greenwich exact - !--------------------------------- - - where(abs(grid_global(:,:,:,1)) < 1.d-10) grid_global(:,:,:,1) = 0.0 - - !--------------------------------- - ! Clean Up Corners - !--------------------------------- - grid_global( 1,1:npts+1,:,2)=grid_global(npts+1,1:npts+1,:,1) - grid_global( 1,1:npts+1,:,3)=grid_global(npts+1:1:-1,npts+1,:,1) - grid_global(1:npts+1,npts+1,:,5)=grid_global(1,npts+1:1:-1,:,1) - grid_global(1:npts+1,npts+1,:,6)=grid_global(1:npts+1,1,:,1) - grid_global(1:npts+1, 1,:,3)=grid_global(1:npts+1,npts+1,:,2) - grid_global(1:npts+1, 1,:,4)=grid_global(npts+1,npts+1:1:-1,:,2) - grid_global(npts+1,1:npts+1,:,6)=grid_global(npts+1:1:-1,1,:,2) - grid_global( 1,1:npts+1,:,4)=grid_global(npts+1,1:npts+1,:,3) - grid_global( 1,1:npts+1,:,5)=grid_global(npts+1:1:-1,npts+1,:,3) - grid_global(npts+1,1:npts+1,:,3)=grid_global(1,1:npts+1,:,4) - grid_global(1:npts+1, 1,:,5)=grid_global(1:npts+1,npts+1,:,4) - grid_global(1:npts+1, 1,:,6)=grid_global(npts+1,npts+1:1:-1,:,4) - grid_global( 1,1:npts+1,:,6)=grid_global(npts+1,1:npts+1,:,5) + where(grid_global(:,:,1,:) < 0.) & + grid_global(:,:,1,:) = grid_global(:,:,1,:) + 2.* PI + + ! Keep Equator and Greenwich exact + !--------------------------------- + + where(abs(grid_global(:,:,:,1)) < 1.d-10) grid_global(:,:,:,1) = 0.0 + + !--------------------------------- + ! Clean Up Corners + !--------------------------------- + grid_global( 1,1:npts+1,:,2)=grid_global(npts+1,1:npts+1,:,1) + grid_global( 1,1:npts+1,:,3)=grid_global(npts+1:1:-1,npts+1,:,1) + grid_global(1:npts+1,npts+1,:,5)=grid_global(1,npts+1:1:-1,:,1) + grid_global(1:npts+1,npts+1,:,6)=grid_global(1:npts+1,1,:,1) + grid_global(1:npts+1, 1,:,3)=grid_global(1:npts+1,npts+1,:,2) + grid_global(1:npts+1, 1,:,4)=grid_global(npts+1,npts+1:1:-1,:,2) + grid_global(npts+1,1:npts+1,:,6)=grid_global(npts+1:1:-1,1,:,2) + grid_global( 1,1:npts+1,:,4)=grid_global(npts+1,1:npts+1,:,3) + grid_global( 1,1:npts+1,:,5)=grid_global(npts+1:1:-1,npts+1,:,3) + grid_global(npts+1,1:npts+1,:,3)=grid_global(1,1:npts+1,:,4) + grid_global(1:npts+1, 1,:,5)=grid_global(1:npts+1,npts+1,:,4) + grid_global(1:npts+1, 1,:,6)=grid_global(npts+1,npts+1:1:-1,:,4) + grid_global( 1,1:npts+1,:,6)=grid_global(npts+1,1:npts+1,:,5) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! ! MAT BMA You cannot refer to FV_Atm here because line 286 says this ! + ! ! MAT BMA You cannot refer to FV_Atm here because line 288 says this ! ! ! is only for FV_Atm not allocated ! ! ! ! !------------------------ ! @@ -347,269 +344,267 @@ function AppGridCreateF(IM_WORLD, JM_WORLD, LM, NX, NY, rc) result(esmfgrid) ! endif ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - else - - if (MAPL_AM_I_ROOT()) write(*,*)'AppGridCreate can not make DP grid by itself' - ASSERT_(.false.) - - end if - - end if - -! Fill lat/lons at cell center locations -! Retrieve the coordinates so we can set them - call ESMF_GridGetCoord(esmfgrid, coordDim=1, localDE=0, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - farrayPtr=lons, rc=status) - VERIFY_(STATUS) - - call ESMF_GridGetCoord(esmfgrid, coordDim=2, localDE=0, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - farrayPtr=lats, rc=status) - VERIFY_(STATUS) - -! Save corners into the esmfgrid -!--------------------------- - allocate(gridCornerLons((size(lons,1)+1)*(size(lons,2)+1)), stat=status) - VERIFY_(STATUS) - allocate(gridCornerLats((size(lats,1)+1)*(size(lats,2)+1)), stat=status) - VERIFY_(STATUS) - - idx = 0 - do jg=jsg,jeg+1 - do ig=isg,ieg+1 - i=ig - j=jg-myTile*npts - idx = idx + 1 - if (AppGridStandAlone) then - gridCornerLons(idx) = grid_global(i, j, 1, myTile+1) - gridCornerLats(idx) = grid_global(i, j, 2, myTile+1) - else - gridCornerLons(idx) = FV_Atm(1)%gridstruct%grid(i,j,1) - gridCornerLats(idx) = FV_Atm(1)%gridstruct%grid(i,j,2) - end if - end do - end do - - call ESMF_AttributeSet(esmfgrid, name='GridCornerLons:', & - itemCount = idx, valueList=gridCornerLons, rc=status) - VERIFY_(STATUS) - call ESMF_AttributeSet(esmfgrid, name='GridCornerLats:', & - itemCount = idx, valueList=gridCornerLats, rc=status) - VERIFY_(STATUS) - - deallocate(gridCornerLats, stat=status) - VERIFY_(STATUS) - deallocate(gridCornerLons, stat=status) - VERIFY_(STATUS) - - do jg=jsg,jeg - do ig=isg,ieg - i=ig - j=jg-myTile*npts - if (AppGridStandAlone) then - call cell_center2(grid_global(i,j, 1:2,myTile+1), grid_global(i+1,j, 1:2,myTile+1), & - grid_global(i,j+1,1:2,myTile+1), grid_global(i+1,j+1,1:2,myTile+1), & - alocs) - else - call cell_center2(dble(FV_Atm(1)%gridstruct%grid(i,j, 1:2)), & - dble(FV_Atm(1)%gridstruct%grid(i+1,j, 1:2)), & - dble(FV_Atm(1)%gridstruct%grid(i,j+1,1:2)), & - dble(FV_Atm(1)%gridstruct%grid(i+1,j+1,1:2)), & - alocs) - end if - i=ig-isg+1 - j=jg-jsg+1 - lons(i,j) = alocs(1) - lats(i,j) = alocs(2) - enddo - enddo - - if (AppGridStandAlone) then - DEALLOCGLOB_(grid_global) - end if - call MAPL_MemUtilsWrite(VM, trim(Iam), RC=STATUS ) - VERIFY_(STATUS) - - RETURN_(ESMF_SUCCESS) + else + + if (MAPL_AM_I_ROOT()) write(*,*)'AppGridCreate can not make DP grid by itself' + _ASSERT(.false.,'needs informative message') + + end if + + end if + + ! Fill lat/lons at cell center locations + ! Retrieve the coordinates so we can set them + call ESMF_GridGetCoord(esmfgrid, coordDim=1, localDE=0, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + farrayPtr=lons, rc=status) + VERIFY_(STATUS) + + call ESMF_GridGetCoord(esmfgrid, coordDim=2, localDE=0, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + farrayPtr=lats, rc=status) + VERIFY_(STATUS) + + ! Save corners into the esmfgrid + !--------------------------- + allocate(gridCornerLons((size(lons,1)+1)*(size(lons,2)+1)), stat=status) + VERIFY_(STATUS) + allocate(gridCornerLats((size(lats,1)+1)*(size(lats,2)+1)), stat=status) + VERIFY_(STATUS) + + idx = 0 + do jg=jsg,jeg+1 + do ig=isg,ieg+1 + i=ig + j=jg-myTile*npts + idx = idx + 1 + if (AppGridStandAlone) then + gridCornerLons(idx) = grid_global(i, j, 1, myTile+1) + gridCornerLats(idx) = grid_global(i, j, 2, myTile+1) + else + gridCornerLons(idx) = FV_Atm(1)%gridstruct%grid(i,j,1) + gridCornerLats(idx) = FV_Atm(1)%gridstruct%grid(i,j,2) + end if + end do + end do + + call ESMF_AttributeSet(esmfgrid, name='GridCornerLons:', & + itemCount = idx, valueList=gridCornerLons, rc=status) + VERIFY_(STATUS) + call ESMF_AttributeSet(esmfgrid, name='GridCornerLats:', & + itemCount = idx, valueList=gridCornerLats, rc=status) + VERIFY_(STATUS) + + deallocate(gridCornerLats, stat=status) + VERIFY_(STATUS) + deallocate(gridCornerLons, stat=status) + VERIFY_(STATUS) + + do jg=jsg,jeg + do ig=isg,ieg + i=ig + j=jg-myTile*npts + if (AppGridStandAlone) then + call cell_center2(grid_global(i,j, 1:2,myTile+1), grid_global(i+1,j, 1:2,myTile+1), & + grid_global(i,j+1,1:2,myTile+1), grid_global(i+1,j+1,1:2,myTile+1), & + alocs) + else + call cell_center2(dble(FV_Atm(1)%gridstruct%grid(i,j, 1:2)), & + dble(FV_Atm(1)%gridstruct%grid(i+1,j, 1:2)), & + dble(FV_Atm(1)%gridstruct%grid(i,j+1,1:2)), & + dble(FV_Atm(1)%gridstruct%grid(i+1,j+1,1:2)), & + alocs) + end if + i=ig-isg+1 + j=jg-jsg+1 + lons(i,j) = alocs(1) + lats(i,j) = alocs(2) + enddo + enddo + + if (AppGridStandAlone) then + DEALLOCGLOB_(grid_global) + end if + call MAPL_MemUtilsWrite(VM, trim(Iam), RC=STATUS ) + VERIFY_(STATUS) + + RETURN_(ESMF_SUCCESS) end function AppGridCreateF subroutine AppGridCreate (META, esmfgrid, RC) #include "MAPL_Generic.h" - use ESMF - use MAPL_Mod - use MAPL_BaseMod - use MAPL_GenericMod - use fv_arrays_mod, only: REAL4, REAL8, R_GRID - implicit none - -! !ARGUMENTS: - - type(MAPL_MetaComp), intent(INOUT) :: META - type (ESMF_Grid), intent( OUT) :: esmfgrid - integer, optional, intent( OUT) :: rc - -! ErrLog variables -!----------------- - - integer :: STATUS - character(len=ESMF_MAXSTR), parameter :: Iam="AppGridCreate" - -! Local variables -!----------------- - - integer :: IM_WORLD - integer :: JM_WORLD - integer :: LM - integer :: NX, NY - character(len=ESMF_MAXSTR) :: tmpname, gridname - character(len=ESMF_MAXSTR) :: FMT, FMTIM, FMTJM - real(REAL8) :: deglon=0.0 - - interface - function AppGridCreateF(IM_WORLD, JM_WORLD, LM, NX, NY, rc) - use ESMF - implicit none - type(ESMF_Grid) :: AppGridCreateF - integer, intent(IN) :: IM_WORLD, JM_WORLD, LM - integer, intent(IN) :: NX, NY - integer, optional, intent(OUT) :: rc - end function AppGridCreateF - end interface - -! Get Decomposition from CF -!-------------------------- - -! !RESOURCE_ITEM: none :: Processing elements in 1st dimension - call MAPL_GetResource( META, NX, label ='NX:', default=1, rc = status ) - VERIFY_(STATUS) -! !RESOURCE_ITEM: none :: Processing elements in 2nd dimension - call MAPL_GetResource( META, NY, label ='NY:', default=1, rc = status ) - VERIFY_(STATUS) - -! Get World problem size from CF -!------------------------------- - -! !RESOURCE_ITEM: none :: Grid size in 1st dimension - call MAPL_GetResource( META, IM_WORLD, 'AGCM_IM:', rc = status ) - VERIFY_(STATUS) -! !RESOURCE_ITEM: none :: Grid size in 2nd dimension - call MAPL_GetResource( META, JM_WORLD, 'AGCM_JM:', rc = status ) - VERIFY_(STATUS) -! !RESOURCE_ITEM: none :: Grid size in 3rd dimension - call MAPL_GetResource( META, LM, 'AGCM_LM:', default=1, rc = status ) - VERIFY_(STATUS) - -! The grid's name is optional -!---------------------------- - - call GET_INT_FORMAT(IM_WORLD, FMTIM) - call GET_INT_FORMAT(JM_WORLD, FMTJM) - FMT = '(A,' // trim(FMTIM) //',A,' // trim(FMTJM) // ',A)' - write(tmpname,trim(FMT)) 'PE',IM_WORLD,'x',JM_WORLD,'-CF' -! !RESOURCE_ITEM: none :: Optional grid name - call MAPL_GetResource( META, GRIDNAME, 'AGCM_GRIDNAME:', default=trim(tmpname), rc = status ) - VERIFY_(STATUS) - - esmfgrid = AppGridCreateF(IM_WORLD, JM_WORLD, LM, NX, NY, STATUS) - VERIFY_(STATUS) - - RETURN_(ESMF_SUCCESS) - - end subroutine AppGridCreate - - subroutine GET_INT_FORMAT(N, FMT) - integer :: N - character(len=*) :: FMT - - IF(N < 10) THEN - FMT = 'I1' - ELSE IF (N< 100) THEN - FMT = 'I2' - ELSE IF (N< 1000) THEN - FMT = 'I3' - ELSE IF (N< 10000) THEN - FMT = 'I4' - else - FMT = 'I5' - end IF - end subroutine GET_INT_FORMAT - - subroutine DecomposeDim ( dim_world,dim,NDEs ) - -! -! From FMS/MPP -! - - implicit none - integer dim_world, NDEs - integer dim(0:NDEs-1) - - integer :: is,ie,isg,ieg - integer :: ndiv,ndivs,imax,ndmax,ndmirror,n - integer :: ibegin(0:NDEs-1) - integer :: iend(0:NDEs-1) - - logical :: symmetrize - logical :: even, odd - even(n) = (mod(n,2).EQ.0) - odd (n) = (mod(n,2).EQ.1) - - isg = 1 - ieg = dim_world - ndivs = NDEs - - is = isg - n = 0 - do ndiv=0,ndivs-1 - !modified for mirror-symmetry - !original line - ! ie = is + CEILING( float(ieg-is+1)/(ndivs-ndiv) ) - 1 - - !problem of dividing nx points into n domains maintaining symmetry - !i.e nx=18 n=4 4554 and 5445 are solutions but 4455 is not. - !this will always work for nx even n even or odd - !this will always work for nx odd, n odd - !this will never work for nx odd, n even: for this case we supersede the mirror calculation - ! symmetrize = .NOT. ( mod(ndivs,2).EQ.0 .AND. mod(ieg-isg+1,2).EQ.1 ) - !nx even n odd fails if n>nx/2 - symmetrize = ( even(ndivs) .AND. even(ieg-isg+1) ) .OR. & - ( odd(ndivs) .AND. odd(ieg-isg+1) ) .OR. & - ( odd(ndivs) .AND. even(ieg-isg+1) .AND. ndivs.LT.(ieg-isg+1)/2 ) - - !mirror domains are stored in the list and retrieved if required. - if( ndiv.EQ.0 )then - !initialize max points and max domains - imax = ieg - ndmax = ndivs - end if - !do bottom half of decomposition, going over the midpoint for odd ndivs - if( ndiv.LT.(ndivs-1)/2+1 )then - !domain is sized by dividing remaining points by remaining domains - ie = is + CEILING( REAL(imax-is+1)/(ndmax-ndiv) ) - 1 - ndmirror = (ndivs-1) - ndiv !mirror domain - if( ndmirror.GT.ndiv .AND. symmetrize )then !only for domains over the midpoint - !mirror extents, the max(,) is to eliminate overlaps - ibegin(ndmirror) = max( isg+ieg-ie, ie+1 ) - iend(ndmirror) = max( isg+ieg-is, ie+1 ) - imax = ibegin(ndmirror) - 1 - ndmax = ndmax - 1 - end if - else - if( symmetrize )then - !do top half of decomposition by retrieving saved values - is = ibegin(ndiv) - ie = iend(ndiv) - else - ie = is + CEILING( REAL(imax-is+1)/(ndmax-ndiv) ) - 1 - end if - end if - dim(ndiv) = ie-is+1 - is = ie + 1 - end do - - end subroutine DecomposeDim + use ESMF + use MAPL + use fv_arrays_mod, only: REAL4, REAL8, R_GRID + implicit none + + ! !ARGUMENTS: + + type(MAPL_MetaComp), intent(INOUT) :: META + type (ESMF_Grid), intent( OUT) :: esmfgrid + integer, optional, intent( OUT) :: rc + + ! ErrLog variables + !----------------- + + integer :: STATUS + character(len=ESMF_MAXSTR), parameter :: Iam="AppGridCreate" + + ! Local variables + !----------------- + + integer :: IM_WORLD + integer :: JM_WORLD + integer :: LM + integer :: NX, NY + character(len=ESMF_MAXSTR) :: tmpname, gridname + character(len=ESMF_MAXSTR) :: FMT, FMTIM, FMTJM + real(REAL8) :: deglon=0.0 + + interface + function AppGridCreateF(IM_WORLD, JM_WORLD, LM, NX, NY, rc) + use ESMF + implicit none + type(ESMF_Grid) :: AppGridCreateF + integer, intent(IN) :: IM_WORLD, JM_WORLD, LM + integer, intent(IN) :: NX, NY + integer, optional, intent(OUT) :: rc + end function AppGridCreateF + end interface + + ! Get Decomposition from CF + !-------------------------- + + ! !RESOURCE_ITEM: none :: Processing elements in 1st dimension + call MAPL_GetResource( META, NX, label ='NX:', default=1, rc = status ) + VERIFY_(STATUS) + ! !RESOURCE_ITEM: none :: Processing elements in 2nd dimension + call MAPL_GetResource( META, NY, label ='NY:', default=1, rc = status ) + VERIFY_(STATUS) + + ! Get World problem size from CF + !------------------------------- + + ! !RESOURCE_ITEM: none :: Grid size in 1st dimension + call MAPL_GetResource( META, IM_WORLD, 'AGCM_IM:', rc = status ) + VERIFY_(STATUS) + ! !RESOURCE_ITEM: none :: Grid size in 2nd dimension + call MAPL_GetResource( META, JM_WORLD, 'AGCM_JM:', rc = status ) + VERIFY_(STATUS) + ! !RESOURCE_ITEM: none :: Grid size in 3rd dimension + call MAPL_GetResource( META, LM, 'AGCM_LM:', default=1, rc = status ) + VERIFY_(STATUS) + + ! The grid's name is optional + !---------------------------- + + call GET_INT_FORMAT(IM_WORLD, FMTIM) + call GET_INT_FORMAT(JM_WORLD, FMTJM) + FMT = '(A,' // trim(FMTIM) //',A,' // trim(FMTJM) // ',A)' + write(tmpname,trim(FMT)) 'PE',IM_WORLD,'x',JM_WORLD,'-CF' + ! !RESOURCE_ITEM: none :: Optional grid name + call MAPL_GetResource( META, GRIDNAME, 'AGCM_GRIDNAME:', default=trim(tmpname), rc = status ) + VERIFY_(STATUS) + + esmfgrid = AppGridCreateF(IM_WORLD, JM_WORLD, LM, NX, NY, STATUS) + VERIFY_(STATUS) + + RETURN_(ESMF_SUCCESS) + +end subroutine AppGridCreate + +subroutine GET_INT_FORMAT(N, FMT) + integer :: N + character(len=*) :: FMT + + IF(N < 10) THEN + FMT = 'I1' + ELSE IF (N< 100) THEN + FMT = 'I2' + ELSE IF (N< 1000) THEN + FMT = 'I3' + ELSE IF (N< 10000) THEN + FMT = 'I4' + else + FMT = 'I5' + end IF +end subroutine GET_INT_FORMAT + +subroutine DecomposeDim ( dim_world,dim,NDEs ) + + ! + ! From FMS/MPP + ! + + implicit none + integer dim_world, NDEs + integer dim(0:NDEs-1) + + integer :: is,ie,isg,ieg + integer :: ndiv,ndivs,imax,ndmax,ndmirror,n + integer :: ibegin(0:NDEs-1) + integer :: iend(0:NDEs-1) + + logical :: symmetrize + logical :: even, odd + even(n) = (mod(n,2).EQ.0) + odd (n) = (mod(n,2).EQ.1) + + isg = 1 + ieg = dim_world + ndivs = NDEs + + is = isg + n = 0 + do ndiv=0,ndivs-1 + !modified for mirror-symmetry + !original line + ! ie = is + CEILING( float(ieg-is+1)/(ndivs-ndiv) ) - 1 + + !problem of dividing nx points into n domains maintaining symmetry + !i.e nx=18 n=4 4554 and 5445 are solutions but 4455 is not. + !this will always work for nx even n even or odd + !this will always work for nx odd, n odd + !this will never work for nx odd, n even: for this case we supersede the mirror calculation + ! symmetrize = .NOT. ( mod(ndivs,2).EQ.0 .AND. mod(ieg-isg+1,2).EQ.1 ) + !nx even n odd fails if n>nx/2 + symmetrize = ( even(ndivs) .AND. even(ieg-isg+1) ) .OR. & + ( odd(ndivs) .AND. odd(ieg-isg+1) ) .OR. & + ( odd(ndivs) .AND. even(ieg-isg+1) .AND. ndivs.LT.(ieg-isg+1)/2 ) + + !mirror domains are stored in the list and retrieved if required. + if( ndiv.EQ.0 )then + !initialize max points and max domains + imax = ieg + ndmax = ndivs + end if + !do bottom half of decomposition, going over the midpoint for odd ndivs + if( ndiv.LT.(ndivs-1)/2+1 )then + !domain is sized by dividing remaining points by remaining domains + ie = is + CEILING( REAL(imax-is+1)/(ndmax-ndiv) ) - 1 + ndmirror = (ndivs-1) - ndiv !mirror domain + if( ndmirror.GT.ndiv .AND. symmetrize )then !only for domains over the midpoint + !mirror extents, the max(,) is to eliminate overlaps + ibegin(ndmirror) = max( isg+ieg-ie, ie+1 ) + iend(ndmirror) = max( isg+ieg-is, ie+1 ) + imax = ibegin(ndmirror) - 1 + ndmax = ndmax - 1 + end if + else + if( symmetrize )then + !do top half of decomposition by retrieving saved values + is = ibegin(ndiv) + ie = iend(ndiv) + else + ie = is + CEILING( REAL(imax-is+1)/(ndmax-ndiv) ) - 1 + end if + end if + dim(ndiv) = ie-is+1 + is = ie + 1 + end do + +end subroutine DecomposeDim diff --git a/CMakeLists.txt b/CMakeLists.txt index d2527e0..fca1d26 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,15 +9,15 @@ set (srcs AdvCore_GridCompMod.F90 DynCore_GridCompMod.F90 CreateInterpWeights_GridCompMod.F90 StandAlone_DynAdvCore_GridCompMod.F90 - CubedSphereGridFactory.F90 CubeToLatLonRegridder.F90 LatLonToCubeRegridder.F90 CubeToCubeRegridder.F90 - CubeToLatLon.F90 -#***DAS***# CubeGridPrototype.F90 + CubeToLatLon.F90 + CubeGridPrototype.F90 + GEOS_FV3_Utilities.F90 ) -set(dependencies MAPL_Base GMAO_hermes GEOS_Shared) +set(dependencies MAPL gftl-shared GMAO_hermes GEOS_Shared) if (FV_PRECISION STREQUAL R4) set (GFDL fms_r4) @@ -63,7 +63,7 @@ set_target_properties (StandAlone_FV3_Dycore.x PROPERTIES LINK_FLAGS "${OpenMP_F ecbuild_add_executable ( TARGET rs_scale.x SOURCES rs_scale.F90 - LIBS ${this} GMAO_pFIO) + LIBS ${this}) set_target_properties (rs_scale.x PROPERTIES LINK_FLAGS "${OpenMP_Fortran_FLAGS}") ecbuild_add_executable ( diff --git a/CreateInterpWeights.F90 b/CreateInterpWeights.F90 index c55e49a..1159a61 100755 --- a/CreateInterpWeights.F90 +++ b/CreateInterpWeights.F90 @@ -4,7 +4,7 @@ program CreateInterpWeights !--------------------------------------------------------------------! ! purpose: driver for MPI-IO test module ! !--------------------------------------------------------------------! - use MAPL_Mod + use MAPL use CreateInterpWeights_GridCompMod, only: SetServices implicit none diff --git a/CreateInterpWeights_GridCompMod.F90 b/CreateInterpWeights_GridCompMod.F90 index 250c001..e90412f 100755 --- a/CreateInterpWeights_GridCompMod.F90 +++ b/CreateInterpWeights_GridCompMod.F90 @@ -13,7 +13,7 @@ Module CreateInterpWeights_GridCompMod ! !USES: use ESMF ! ESMF base class - use MAPL_Mod ! GEOS base class + use MAPL ! GEOS base class ! !PUBLIC MEMBER FUNCTIONS: diff --git a/Cube2LatLon.F90 b/Cube2LatLon.F90 index 54a8439..73ab820 100644 --- a/Cube2LatLon.F90 +++ b/Cube2LatLon.F90 @@ -3,7 +3,7 @@ subroutine cube2latlon(npx, npy, nlon, nlat, data_cs, data_ll) #define REAL8 8 use ESMF - use MAPL_Mod, only : MAPL_UNDEF + use MAPL, only : MAPL_UNDEF use MAPL_IOMod, only : GETFILEUNIT, FREE_FILE use CUB2LATLON_mod, only : init_latlon_grid, & read_c2l_weight, write_c2l_weight, & diff --git a/CubeGridPrototype.F90 b/CubeGridPrototype.F90 new file mode 100644 index 0000000..be70886 --- /dev/null +++ b/CubeGridPrototype.F90 @@ -0,0 +1,27 @@ +module CubeGridPrototype + implicit none + + private + public :: register_grid_and_regridders +contains + + subroutine register_grid_and_regridders() + use MAPL_RegridderManagerMod, only: regridder_manager + use MAPL_RegridderSpecMod, only: REGRID_METHOD_BILINEAR + use LatLonToCubeRegridderMod + use CubeToLatLonRegridderMod + use CubeToCubeRegridderMod + + type (CubeToLatLonRegridder) :: cube_to_latlon_prototype + type (LatLonToCubeRegridder) :: latlon_to_cube_prototype + type (CubeToCubeRegridder) :: cube_to_cube_prototype + + associate (method => REGRID_METHOD_BILINEAR, mgr => regridder_manager) + call mgr%add_prototype('Cubed-Sphere', 'LatLon', method, cube_to_latlon_prototype) + call mgr%add_prototype('LatLon', 'Cubed-Sphere', method, latlon_to_cube_prototype) + call mgr%add_prototype('Cubed-Sphere', 'Cubed-Sphere', method, cube_to_cube_prototype) + end associate + + end subroutine register_grid_and_regridders + +end module CubeGridPrototype diff --git a/CubeToCubeRegridder.F90 b/CubeToCubeRegridder.F90 index a8eb0c3..8b5c619 100644 --- a/CubeToCubeRegridder.F90 +++ b/CubeToCubeRegridder.F90 @@ -5,10 +5,8 @@ #define _RETURN(A) if(present(rc)) rc=A; return module CubeToCubeRegridderMod - use MAPL_AbstractRegridderMod + use MAPL use CubeLatLonTransformMod - use MAPL_GridSpecMod - use MAPL_RegridderSpecMod use, intrinsic :: iso_fortran_env, only: REAL32 use, intrinsic :: iso_fortran_env, only: REAL64 implicit none @@ -47,9 +45,7 @@ end function newCubeToCubeRegridder subroutine initialize_subclass(this, unusable, rc) - use MAPL_KeywordEnforcerMod - use MAPL_BaseMod - use MAPL_CommsMod + use MAPL use ESMF class (CubeToCubeRegridder), intent(inout) :: this class (KeywordEnforcer), optional, intent(in) :: unusable @@ -115,8 +111,7 @@ end subroutine regrid_scalar_2d_real32 subroutine regrid_scalar_3d_real32(this, q_in, q_out, rc) - use MAPL_CommsMod - use MAPL_BaseMod + use MAPL class (CubeToCubeRegridder), intent(in) :: this real (kind=REAL32), intent(in) :: q_in(:,:,:) diff --git a/CubeToLatLon.F90 b/CubeToLatLon.F90 index edc8b33..913ad95 100644 --- a/CubeToLatLon.F90 +++ b/CubeToLatLon.F90 @@ -1,17 +1,14 @@ ! $Id$ -#define SUCCESS 0 -#define VERIFY_(A) if((A)/=0) then; if(present(rc)) rc=A; PRINT *, Iam, __LINE__; return; endif -#define ASSERT_(A) if(.not.(A)) then; if(present(rc)) rc=1; PRINT *, Iam, __LINE__; return; endif -#define RETURN_(A) if(present(rc)) rc=A; return +#include "MAPL_ErrLog.h" #define DEALLOCGLOB_(A) call deallocGlob(A,status);VERIFY_(status) #define DEALLOCLOCL_(A) if(associated(A)) then; deallocate(A, stat=STATUS); VERIFY_(STATUS); NULLIFY(A); endif #ifdef TAU_PROFILE -#undef ASSERT_ -#define ASSERT_(A) +#undef _ASSERT +#define _ASSERT(A,'needs informative message') #undef VERIFY_ #define VERIFY_(A) @@ -23,11 +20,7 @@ Module CubeLatLonTransformMod use ESMF - use MAPL_BaseMod - use MAPL_LocStreamMod - use MAPL_IOMod - use MAPL_CommsMod - use MAPL_ShmemMod + use MAPL use, intrinsic :: iso_fortran_env, only: REAL64 implicit none @@ -251,7 +244,7 @@ subroutine CubeLatLonDestroy( Tr, rc) ! unfortunately MAPL does not destroy LocationStreams yet end if - RETURN_(SUCCESS) + RETURN_(_SUCCESS) end subroutine CubeLatLonDestroy logical function CubeLatLonIsCreated(Tr) @@ -313,7 +306,7 @@ function CubeLatLonCreate( npx, npy, nlon, nlat, lons, lats, local_ij, rc ) resu ! Begin !------ - ASSERT_(.not.Tr%Created) + _ASSERT(.not.Tr%Created,'needs informative message') npts = npx + 1 @@ -344,28 +337,46 @@ function CubeLatLonCreate( npx, npy, nlon, nlat, lons, lats, local_ij, rc ) resu DEALLOCLOCL_(Tr%g1) DEALLOCLOCL_(Tr%g2) - call MAPL_AllocNodeArray(Tr%index,(/3,nlon,nlat/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%index(3,nlon,nlat),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%index,(/3,nlon,nlat/),rc=STATUS) + else + allocate(Tr%index(3,nlon,nlat),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%weight,(/4,nlon,nlat/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%weight(4,nlon,nlat),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%weight,(/4,nlon,nlat/),rc=STATUS) + else + allocate(Tr%weight(4,nlon,nlat),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%l2c,(/4,npx,npy/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%l2c(4,npx,npy),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%l2c,(/4,npx,npy/),rc=STATUS) + else + allocate(Tr%l2c(4,npx,npy),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%id1,(/npx,npy/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%id1(npx,npy),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%id1,(/npx,npy/),rc=STATUS) + else + allocate(Tr%id1(npx,npy),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%id2,(/npx,npy/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%id2(npx,npy),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%id2,(/npx,npy/),rc=STATUS) + else + allocate(Tr%id2(npx,npy),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%jdc,(/npx,npy/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%jdc(npx,npy),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%jdc,(/npx,npy/),rc=STATUS) + else + allocate(Tr%jdc(npx,npy),stat=status) + end if VERIFY_(STATUS) allocate(Tr%elon(size(lons),size(lats),3),stat=STATUS) @@ -389,28 +400,46 @@ function CubeLatLonCreate( npx, npy, nlon, nlat, lons, lats, local_ij, rc ) resu VERIFY_(STATUS) tr%elat_local => tr%elat(i0:i1,j0:j1,:) - call MAPL_AllocNodeArray(Tr%ee1,(/npx,npy,3/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%ee1(npx,npy,3),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%ee1,(/npx,npy,3/),rc=STATUS) + else + allocate(Tr%ee1(npx,npy,3),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%ee2,(/npx,npy,3/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%ee2(npx,npy,3),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%ee2,(/npx,npy,3/),rc=STATUS) + else + allocate(Tr%ee2(npx,npy,3),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%ff1,(/npx,npy,3/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%ff1(npx,npy,3),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%ff1,(/npx,npy,3/),rc=STATUS) + else + allocate(Tr%ff1(npx,npy,3),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%ff2,(/npx,npy,3/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%ff2(npx,npy,3),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%ff2,(/npx,npy,3/),rc=STATUS) + else + allocate(Tr%ff2(npx,npy,3),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%gg1,(/npx,npy,3/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%gg1(npx,npy,3),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%gg1,(/npx,npy,3/),rc=STATUS) + else + allocate(Tr%gg1(npx,npy,3),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%gg2,(/npx,npy,3/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%gg2(npx,npy,3),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%gg2,(/npx,npy,3/),rc=STATUS) + else + allocate(Tr%gg2(npx,npy,3),stat=status) + end if VERIFY_(STATUS) ! Argument AmNodeRoot passed to GetWeights identifies if we're using SHMEM @@ -468,7 +497,7 @@ function CubeLatLonCreate( npx, npy, nlon, nlat, lons, lats, local_ij, rc ) resu Tr%Created=.true. - RETURN_(SUCCESS) + RETURN_(_SUCCESS) end function CubeLatLonCreate subroutine CubeCubeDestroy( Tr, rc) @@ -486,7 +515,7 @@ subroutine CubeCubeDestroy( Tr, rc) Tr%Created = .false. - RETURN_(SUCCESS) + RETURN_(_SUCCESS) end subroutine CubeCubeDestroy logical function CubeCubeIsCreated(Tr) @@ -516,7 +545,7 @@ function CubeCubeCreate( npx, npy, npxout, npyout, rc ) result(Tr) ! Begin !------ - ASSERT_(.not.Tr%Created) + _ASSERT(.not.Tr%Created,'needs informative message') !ALT npts = npx + 1 npts = npxout ! + 1 @@ -539,30 +568,48 @@ function CubeCubeCreate( npx, npy, npxout, npyout, rc ) result(Tr) DEALLOCGLOB_(Tr%ff2) ! ALT: index and weight are allocated at the output grid resolution - call MAPL_AllocNodeArray(Tr%weight,(/4,npxout,npyout/6,6/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%weight(4,npxout,npyout/6,6),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%weight,(/4,npxout,npyout/6,6/),rc=STATUS) + else + allocate(Tr%weight(4,npxout,npyout/6,6),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%index,(/3,npxout,npyout/6,6/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%index(3,npxout,npyout/6,6),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%index,(/3,npxout,npyout/6,6/),rc=STATUS) + else + allocate(Tr%index(3,npxout,npyout/6,6),stat=status) + end if VERIFY_(STATUS) ! ALT: ff1 and ff2 are allocated at the input grid resolution - call MAPL_AllocNodeArray(Tr%ff1,(/npx,npy,3/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%ff1(npx,npy,3),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%ff1,(/npx,npy,3/),rc=STATUS) + else + allocate(Tr%ff1(npx,npy,3),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%ff2,(/npx,npy,3/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%ff2(npx,npy,3),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%ff2,(/npx,npy,3/),rc=STATUS) + else + allocate(Tr%ff2(npx,npy,3),stat=status) + end if VERIFY_(STATUS) ! ALT: ee1 and ee2 are allocated at the output grid resolution - call MAPL_AllocNodeArray(Tr%ee1,(/npxout,npyout,3/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%ee1(npxout,npyout,3),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%ee1,(/npxout,npyout,3/),rc=STATUS) + else + allocate(Tr%ee1(npxout,npyout,3),stat=status) + end if VERIFY_(STATUS) - call MAPL_AllocNodeArray(Tr%ee2,(/npxout,npyout,3/),rc=STATUS) - if(STATUS==MAPL_NoShm) allocate(Tr%ee2(npxout,npyout,3),stat=status) + if(MAPL_ShmInitialized) then + call MAPL_AllocNodeArray(Tr%ee2,(/npxout,npyout,3/),rc=STATUS) + else + allocate(Tr%ee2(npxout,npyout,3),stat=status) + end if VERIFY_(STATUS) if (MAPL_AmNodeRoot .or. (.not. MAPL_ShmInitialized)) then @@ -575,7 +622,7 @@ function CubeCubeCreate( npx, npy, npxout, npyout, rc ) result(Tr) Tr%Created=.true. - RETURN_(SUCCESS) + RETURN_(_SUCCESS) end function CubeCubeCreate subroutine CubeToCube(Tr, data_cs_in, data_cs_out, rc) @@ -591,7 +638,7 @@ subroutine CubeToCube(Tr, data_cs_in, data_cs_out, rc) integer :: nx,j1,j2,status,itile real(REAL64), allocatable :: var_cs_in(:,:,:), var_cs_out(:,:,:) - ASSERT_(Tr%Created) + _ASSERT(Tr%Created,'needs informative message') nx = Tr%npx @@ -626,7 +673,7 @@ subroutine CubeToCube(Tr, data_cs_in, data_cs_out, rc) deallocate ( var_cs_in, var_cs_out, stat=status) VERIFY_(STATUS) - RETURN_(SUCCESS) + RETURN_(_SUCCESS) end subroutine CubeToCube subroutine C2CInterp(var_in, var_out, index_c2c, weight_c2c) @@ -759,7 +806,7 @@ subroutine CubeToLatLonr8( Tr, data_cs, data_ll, transpose, misval, rc) integer :: nx,j1,j2,status,itile real(REAL64), allocatable :: var_cs(:,:,:) - ASSERT_(Tr%Created) + _ASSERT(Tr%Created,'needs informative message') nx = Tr%npx @@ -770,7 +817,7 @@ subroutine CubeToLatLonr8( Tr, data_cs, data_ll, transpose, misval, rc) allocate ( var_cs(0:nx+1,0:nx+1,ntiles),stat=status) VERIFY_(STATUS) - ASSERT_(.not. (transpose .and. present(misval))) + _ASSERT(.not. (transpose .and. present(misval)),'needs informative message') var_cs=0.0 @@ -797,7 +844,7 @@ subroutine CubeToLatLonr8( Tr, data_cs, data_ll, transpose, misval, rc) deallocate ( var_cs ,stat=status) VERIFY_(STATUS) - RETURN_(SUCCESS) + RETURN_(_SUCCESS) end subroutine CubeToLatLonr8 subroutine CubeToLatLonr4( Tr, data_cs, data_ll, transpose, misval, rc) @@ -819,7 +866,7 @@ subroutine CubeToLatLonr4( Tr, data_cs, data_ll, transpose, misval, rc) real(REAL64), allocatable :: data_cs8(:,:) - ASSERT_(Tr%Created) + _ASSERT(Tr%Created,'needs informative message') nx = Tr%npx @@ -837,7 +884,7 @@ subroutine CubeToLatLonr4( Tr, data_cs, data_ll, transpose, misval, rc) VERIFY_(STATUS) endif - ASSERT_(.not. (transpose .and. present(misval))) + _ASSERT(.not. (transpose .and. present(misval)),'needs informative message') var_cs=0.0 @@ -884,7 +931,7 @@ subroutine CubeToLatLonr4( Tr, data_cs, data_ll, transpose, misval, rc) deallocate ( var_cs, data_ll8 ) - RETURN_(SUCCESS) + RETURN_(_SUCCESS) end subroutine CubeToLatLonr4 subroutine LatLonToCuber8( Tr, data_ll, data_cs, transpose, misval, rc) @@ -902,7 +949,7 @@ subroutine LatLonToCuber8( Tr, data_ll, data_cs, transpose, misval, rc) integer :: nx,j1,j2,status,itile real(REAL64), allocatable :: var_cs(:,:,:) - ASSERT_(Tr%Created) + _ASSERT(Tr%Created,'needs informative message') nx = Tr%npx @@ -938,7 +985,7 @@ subroutine LatLonToCuber8( Tr, data_ll, data_cs, transpose, misval, rc) deallocate ( var_cs ,STAT=STATUS) VERIFY_(STATUS) - RETURN_(SUCCESS) + RETURN_(_SUCCESS) end subroutine LatLonToCuber8 subroutine LatLonToCuber4( Tr, data_ll, data_cs, transpose, misval, rc) @@ -960,7 +1007,7 @@ subroutine LatLonToCuber4( Tr, data_ll, data_cs, transpose, misval, rc) !JK for conservative interp-------------- - ASSERT_(Tr%Created) + _ASSERT(Tr%Created,'needs informative message') nx = Tr%npx @@ -1033,7 +1080,7 @@ subroutine LatLonToCuber4( Tr, data_ll, data_cs, transpose, misval, rc) deallocate ( data_cs8 ) deallocate ( data_ll8 ) - RETURN_(SUCCESS) + RETURN_(_SUCCESS) end subroutine LatLonToCuber4 subroutine C2LInterp(cubsph, latlon, index, weight, misval, subset, transpose) @@ -1364,7 +1411,7 @@ subroutine SphericalToCartesianR4(Tr, U, V, Uxyz, Transpose, SphIsLL, Rotate, RC real(REAL64), pointer :: e1(:,:,:), e2(:,:,:) if(.not.Rotate) then - ASSERT_(.not.Transpose .and. .not.SphIsLL) + _ASSERT(.not.Transpose .and. .not.SphIsLL,'needs informative message') end if if(SphIsLL) then @@ -1437,7 +1484,7 @@ subroutine SphericalToCartesianREAL64(Tr, U, V, Uxyz, Transpose, SphIsLL, Rotate real(REAL64), pointer :: e1(:,:,:), e2(:,:,:) if(.not.Rotate) then - ASSERT_(.not.Transpose .and. .not.SphIsLL) + _ASSERT(.not.Transpose .and. .not.SphIsLL,'needs informative message') end if if(SphIsLL) then @@ -1510,7 +1557,7 @@ subroutine CartesianToSphericalR4(Tr, Uxyz, U, V, Transpose, SphIsLL, Rotate, RC real(REAL64), pointer :: e1(:,:,:), e2(:,:,:) if(.not.Rotate) then - ASSERT_(.not.Transpose .and. .not.SphIsLL) + _ASSERT(.not.Transpose .and. .not.SphIsLL,'needs informative message') end if if(SphIsLL) then @@ -1588,7 +1635,7 @@ subroutine CartesianToSphericalREAL64(Tr, Uxyz, U, V, Transpose, SphIsLL, Rotate real(REAL64), pointer :: e1(:,:,:), e2(:,:,:) if(.not.Rotate) then - ASSERT_(.not.Transpose .and. .not.SphIsLL) + _ASSERT(.not.Transpose .and. .not.SphIsLL,'needs informative message') end if if(SphIsLL) then @@ -2094,8 +2141,11 @@ subroutine deallocGlob_i4_4(a, status) if (associated(a)) then a=0 - call MAPL_DeAllocNodeArray(a, rc=status) - if (status==MAPL_NoShm) deallocate(a, stat=status) + if(MAPL_ShmInitialized) then + call MAPL_DeAllocNodeArray(a, rc=status) + else + deallocate(a, stat=status) + end if if (status /= 0) return nullify(a) else @@ -2110,8 +2160,11 @@ subroutine deallocGlob_i4_3(a, status) if (associated(a)) then a=0 - call MAPL_DeAllocNodeArray(a, rc=status) - if (status==MAPL_NoShm) deallocate(a, stat=status) + if(MAPL_ShmInitialized) then + call MAPL_DeAllocNodeArray(a, rc=status) + else + deallocate(a, stat=status) + end if if (status /= 0) return nullify(a) else @@ -2126,8 +2179,11 @@ subroutine deallocGlob_i4_2(a, status) if (associated(a)) then a=0 - call MAPL_DeAllocNodeArray(a, rc=status) - if (status==MAPL_NoShm) deallocate(a, stat=status) + if(MAPL_ShmInitialized) then + call MAPL_DeAllocNodeArray(a, rc=status) + else + deallocate(a, stat=status) + end if if (status /= 0) return nullify(a) else @@ -2142,8 +2198,11 @@ subroutine deallocGlob_r8_3(a, status) if (associated(a)) then a=0 - call MAPL_DeAllocNodeArray(a, rc=status) - if (status==MAPL_NoShm) deallocate(a, stat=status) + if(MAPL_ShmInitialized) then + call MAPL_DeAllocNodeArray(a, rc=status) + else + deallocate(a, stat=status) + end if if (status /= 0) return nullify(a) else @@ -2158,8 +2217,11 @@ subroutine deallocGlob_r8_4(a, status) if (associated(a)) then a=0 - call MAPL_DeAllocNodeArray(a, rc=status) - if (status==MAPL_NoShm) deallocate(a, stat=status) + if(MAPL_ShmInitialized) then + call MAPL_DeAllocNodeArray(a, rc=status) + else + deallocate(a, stat=status) + end if if (status /= 0) return nullify(a) else diff --git a/CubeToLatLonRegridder.F90 b/CubeToLatLonRegridder.F90 index 8c64602..3b552b9 100644 --- a/CubeToLatLonRegridder.F90 +++ b/CubeToLatLonRegridder.F90 @@ -5,10 +5,8 @@ #define _RETURN(A) if(present(rc)) rc=A; return module CubeToLatLonRegridderMod - use MAPL_AbstractRegridderMod + use MAPL use CubeLatLonTransformMod - use MAPL_GridSpecMod - use MAPL_RegridderSpecMod use, intrinsic :: iso_fortran_env, only: REAL32 use, intrinsic :: iso_fortran_env, only: REAL64 implicit none @@ -47,11 +45,7 @@ function newCubeToLatLonRegridder(regridder_spec, rc) result(regridder) end function newCubeToLatLonRegridder subroutine initialize_subclass(this, unusable, rc) - use MAPL_KeywordEnforcerMod - use MAPL_BaseMod - use MAPL_AbstractGridFactoryMod - use MAPL_LatLonGridFactoryMod - use MAPL_GridManagerMod + use MAPL class (CubeToLatLonRegridder), intent(inout) :: this class (KeywordEnforcer), optional, intent(in) :: unusable integer, optional, intent(out) :: rc @@ -135,8 +129,7 @@ end subroutine regrid_scalar_2d_real32 subroutine regrid_scalar_3d_real32(this, q_in, q_out, rc) - use MAPL_CommsMod - use MAPL_BaseMod + use MAPL class (CubeToLatLonRegridder), intent(in) :: this real (kind=REAL32), intent(in) :: q_in(:,:,:) @@ -218,8 +211,7 @@ end subroutine regrid_scalar_3d_real32 subroutine regrid_vector_3d_real32(this, u_in, v_in, u_out, v_out, rotate, rc) - use MAPL_CommsMod - use MAPL_BaseMod + use MAPL use, intrinsic :: iso_fortran_env, only: REAL32 class (CubeToLatLonRegridder), intent(in) :: this real(kind=REAL32), intent(in) :: u_in(:,:,:) @@ -304,8 +296,7 @@ subroutine transpose_regrid_scalar_2d_real32(this, q_in, q_out, rc) end subroutine transpose_regrid_scalar_2d_real32 subroutine transpose_regrid_scalar_3d_real32(this, q_in, q_out, rc) - use MAPL_CommsMod - use MAPL_BaseMod + use MAPL class (CubeToLatLonRegridder), intent(in) :: this real (kind=REAL32), intent(in) :: q_in(:,:,:) @@ -386,8 +377,7 @@ subroutine transpose_regrid_scalar_3d_real32(this, q_in, q_out, rc) end subroutine transpose_regrid_scalar_3d_real32 subroutine transpose_regrid_vector_3d_real32(this, u_in, v_in, u_out, v_out, rotate, rc) - use MAPL_CommsMod - use MAPL_BaseMod + use MAPL use, intrinsic :: iso_fortran_env, only: REAL32 class (CubeToLatLonRegridder), intent(in) :: this real(kind=REAL32), intent(in) :: u_in(:,:,:) diff --git a/CubedSphereGridFactory.F90 b/CubedSphereGridFactory.F90 deleted file mode 100644 index b1b3021..0000000 --- a/CubedSphereGridFactory.F90 +++ /dev/null @@ -1,706 +0,0 @@ -!----------------------------------------------------- -! Note that this implementation only supports -! "square" faces on the cube. I.e. the number of -! cells along each axis (of each face) are the same. -! IM_WORLD is used for this quantity, and there is no -! equivalent for the "other" axis. -!----------------------------------------------------- - -#define _SUCCESS 0 -#define _FAILURE 1 -#define _VERIFY(A) if( A/=0) then; if(present(rc)) rc=A; PRINT *, Iam, __LINE__; return; endif -#define _ASSERT(A) if(.not.(A)) then; if(present(rc)) rc=_FAILURE; PRINT *, Iam, __LINE__; return; endif -#define _RETURN(A) if(present(rc)) rc=A; return -#include "unused_dummy.H" - - -module CubedSphereGridFactoryMod - use fv_grid_utils_mod, only: gnomonic_grids, cell_center2 - use fv_grid_tools_mod, only: mirror_grid - use MAPL_AbstractGridFactoryMod - use MAPL_MinMaxMod - use MAPL_KeywordEnforcerMod - use ESMF - use MAPL_CommsMod - use MAPL_IOMod, only : GETFILE, FREE_FILE - use fms_mod, only: fms_init - use mpp_domains_mod, only: domain2d, mpp_update_domains - use mpp_mod, only: mpp_error, FATAL, NOTE - !use CubedSphereA2D2CMod - use, intrinsic :: iso_fortran_env, only: REAL64 - implicit none - private - - public :: CubedSphereGridFactory - - integer, parameter :: ndims = 2 - integer, parameter :: UNDEFINED_INTEGER = 1-huge(1) - real, parameter :: UNDEFINED_REAL = huge(1.) - character(len=*), parameter :: UNDEFINED_CHAR = '**' - - integer, parameter :: FV_GRID_TYPE_DEFAULT = 0 - character(len=*), parameter :: GRID_NAME_DEFAULT = 'UNKNOWN' - - integer, parameter :: NUM_CUBE_FACES = 6 - - interface - subroutine A2D2C(U,V,npz,getC) - implicit none - real, intent(INOUT) :: U(:,:,:) - real, intent(INOUT) :: V(:,:,:) - integer, intent( IN) :: npz - logical, optional, intent( IN) :: getC - end subroutine A2D2C - end interface - - type, extends(AbstractGridFactory) :: CubedSphereGridFactory - private - - - character(len=:), allocatable :: grid_name - integer :: grid_type = UNDEFINED_INTEGER - - ! Grid dimensions - Note that we only support "square" grids - integer :: im_world = UNDEFINED_INTEGER - integer :: lm = UNDEFINED_INTEGER - integer :: ntiles = NUM_CUBE_FACES - - ! Domain decomposition: - note that we only support "square" dec - integer :: nx = UNDEFINED_INTEGER - integer :: ny = UNDEFINED_INTEGER - integer, allocatable :: ims(:) - integer, allocatable :: jms(:) - ! rectangle decomposition - integer, allocatable :: jms_2d(:,:) - - ! For halo - type (domain2d) :: domain - - logical :: halo_initialized = .false. - - contains - - procedure :: make_new_grid - procedure :: create_basic_grid - - procedure :: initialize_from_config_with_prefix - procedure :: initialize_from_esmf_distGrid - - procedure :: halo_init - procedure :: halo - - procedure :: check_and_fill_consistency - procedure :: equals - procedure :: generate_grid_name - procedure :: to_string - - procedure :: a2d2c => a2d2c_factory - end type CubedSphereGridFactory - - character(len=*), parameter :: MOD_NAME = 'CubedSphereGridFactory::' - - interface CubedSphereGridFactory - module procedure CubedSphereGridFactory_from_parameters - end interface CubedSphereGridFactory - - interface set_with_default - module procedure set_with_default_integer - module procedure set_with_default_real - module procedure set_with_default_character - module procedure set_with_default_bounds - end interface set_with_default - - -contains - - - function CubedSphereGridFactory_from_parameters(unusable, grid_name, grid_type, & - & im_world, lm, nx, ny, ims, jms, & - & rc) result(factory) - type (CubedSphereGridFactory) :: factory - class (KeywordEnforcer), optional, intent(in) :: unusable - character(len=*), optional, intent(in) :: grid_name - integer, optional, intent(in) :: grid_type - - ! grid details: - integer, optional, intent(in) :: im_world - integer, optional, intent(in) :: lm - - ! decomposition: - integer, optional, intent(in) :: nx - integer, optional, intent(in) :: ny - integer, optional, intent(in) :: ims(:) - integer, optional, intent(in) :: jms(:) - - integer, optional, intent(out) :: rc - - integer :: status - character(len=*), parameter :: Iam = MOD_NAME // 'CubedSphereGridFactory_from_parameters' - - if (present(unusable)) print*,shape(unusable) - - call set_with_default(factory%grid_name, grid_name, GRID_NAME_DEFAULT) - call set_with_default(factory%grid_type, grid_type, FV_GRID_TYPE_DEFAULT) - - call set_with_default(factory%nx, nx, UNDEFINED_INTEGER) - call set_with_default(factory%ny, ny, UNDEFINED_INTEGER) - - call set_with_default(factory%im_world, im_world, UNDEFINED_INTEGER) - call set_with_default(factory%lm, lm, UNDEFINED_INTEGER) - - ! default is unallocated - if (present(ims)) factory%ims = ims - if (present(jms)) factory%jms = jms - - call factory%check_and_fill_consistency(rc=status) - - _VERIFY(status) - - _RETURN(_SUCCESS) - - end function CubedSphereGridFactory_from_parameters - - - function make_new_grid(this, unusable, rc) result(grid) - type (ESMF_Grid) :: grid - class (CubedSphereGridFactory), intent(in) :: this - class (KeywordEnforcer), optional, intent(in) :: unusable - integer, optional, intent(out) :: rc - - integer :: status - character(len=*), parameter :: Iam = MOD_NAME // 'make_grid' - - grid = this%create_basic_grid(rc=status) - _VERIFY(status) - - _RETURN(_SUCCESS) - - end function make_new_grid - - - - function create_basic_grid(this, unusable, rc) result(grid) - type (ESMF_Grid) :: grid - class (CubedSphereGridFactory), intent(in) :: this - class (KeywordEnforcer), optional, intent(in) :: unusable - integer, optional, intent(out) :: rc - - integer :: i,nTile - integer, allocatable :: ims(:,:), jms(:,:) - real(kind=ESMF_KIND_R8), pointer :: centers(:,:) - integer :: status - character(len=*), parameter :: Iam = MOD_NAME // 'create_basic_grid' - - if (this%grid_type <=3) then - nTile=6 - else - nTile=1 - end if - - allocate(ims(this%nx,nTile)) - do i=1,nTile - ims(:,i)=this%ims - enddo - - if(allocated(this%jms_2d)) then - _ASSERT(size(this%jms_2d,2) == 6) - allocate(jms, source = this%jms_2d) - else - allocate(jms(this%ny,nTile)) - do i=1,nTile - jms(:,i)=this%jms - end do - endif - - if (this%grid_type <= 3) then - grid = ESMF_GridCreateCubedSPhere(this%im_world,countsPerDEDim1PTile=ims, & - countsPerDEDim2PTile=jms ,name=this%grid_name, & - staggerLocList=[ESMF_STAGGERLOC_CENTER,ESMF_STAGGERLOC_CORNER], coordSys=ESMF_COORDSYS_SPH_RAD, rc=status) - _VERIFY(status) - call ESMF_AttributeSet(grid, 'GridType', 'Cubed-Sphere', rc=status) - else - grid = ESMF_GridCreateNoPeriDim( & - & name = this%grid_name, & - & countsPerDEDim1=this%ims, & - & countsPerDEDim2=this%jms, & - & indexFlag=ESMF_INDEX_DELOCAL, & - & gridEdgeLWidth=[0,0], & - & gridEdgeUWidth=[1,1], & - & coordDep1=[1,2], & - & coordDep2=[1,2], & - & coordSys=ESMF_COORDSYS_SPH_RAD, & - & rc=status) - _VERIFY(status) - call ESMF_AttributeSet(grid, 'GridType', 'Doubly-Periodic', rc=status) - _VERIFY(status) - call ESMF_GridAddCoord(grid,rc=status) - _VERIFY(status) - call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - farrayPtr=centers, rc=status) - _VERIFY(status) - centers=0.0 - call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - farrayPtr=centers, rc=status) - _VERIFY(status) - centers=0.0 - - end if - - deallocate(ims,jms) - - if (this%lm /= UNDEFINED_INTEGER) then - call ESMF_AttributeSet(grid, name='GRID_LM', value=this%lm, rc=status) - _VERIFY(status) - end if - - call ESMF_AttributeSet(grid, name='NEW_CUBE', value=1,rc=status) - _VERIFY(status) - - _RETURN(_SUCCESS) - end function create_basic_grid - - subroutine initialize_from_config_with_prefix(this, config, prefix, unusable, rc) - use esmf - class (CubedSphereGridFactory), intent(inout) :: this - type (ESMF_Config), intent(inout) :: config - character(len=*), intent(in) :: prefix ! effectively optional due to overload without this argument - class (KeywordEnforcer), optional, intent(in) :: unusable - integer, optional, intent(out) :: rc - - integer :: status - character(len=*), parameter :: Iam = MOD_NAME//'make_geos_grid_from_config' - character(len=ESMF_MAXSTR) :: tmp - type (ESMF_VM) :: vm - integer :: vmcomm, ndes - - if (present(unusable)) print*,shape(unusable) - - call ESMF_VMGetCurrent(vm, rc=status) - _VERIFY(status) - - call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'GRIDNAME:', default=GRID_NAME_DEFAULT) - this%grid_name = trim(tmp) - - call ESMF_ConfigGetAttribute(config, this%grid_type, label=prefix//'CS_GRID_TYPE:', default=FV_GRID_TYPE_DEFAULT) - - call ESMF_ConfigGetAttribute(config, this%nx, label=prefix//'NX:', default=UNDEFINED_INTEGER) - call ESMF_ConfigGetAttribute(config, this%ny, label=prefix//'NY:', default=UNDEFINED_INTEGER) - - call ESMF_ConfigGetAttribute(config, this%im_world, label=prefix//'IM_WORLD:', default=UNDEFINED_INTEGER) - - call get_multi_integer(this%ims, 'IMS:', rc=status) - _VERIFY(status) - - call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'JMS_FILE:', rc=status) - if (status == _SUCCESS) then - call get_jms_from_file(this%jms_2d, trim(tmp),this%ny, rc=status) - _VERIFY(status) - else - call get_multi_integer(this%jms, 'JMS:', rc=status) - _VERIFY(status) - endif - - call ESMF_ConfigGetAttribute(config, this%lm, label=prefix//'LM:', default=UNDEFINED_INTEGER) - - call this%check_and_fill_consistency(rc=status) - _VERIFY(status) - - ! halo initialization - - call ESMF_VmGet(VM, mpicommunicator=vmcomm, petCount=ndes, rc=status) - _VERIFY(status) - - - _RETURN(_SUCCESS) - - contains - - subroutine get_multi_integer(values, label, rc) - integer, allocatable, intent(out) :: values(:) - character(len=*) :: label - integer, optional, intent(out) :: rc - - integer :: i - integer :: n - integer :: tmp - integer :: status - - call ESMF_ConfigFindLabel(config, label=prefix//label,rc=status) - if (status /= _SUCCESS) then - _RETURN(_SUCCESS) - end if - - ! First pass: count values - n = 0 - do - call ESMF_ConfigGetAttribute(config, tmp, default=UNDEFINED_INTEGER, rc=status) - if (status /= _SUCCESS) then - exit - else - n = n + 1 - end if - end do - - ! Second pass: allocate and fill - allocate(values(n), stat=status) ! no point in checking status - _VERIFY(status) - call ESMF_ConfigFindLabel(config, label=prefix//label,rc=status) - do i = 1, n - call ESMF_ConfigGetAttribute(config, values(i), rc=status) - _VERIFY(status) - end do - - _RETURN(_SUCCESS) - - end subroutine get_multi_integer - - subroutine get_jms_from_file(values, file_name, n, rc) - integer, allocatable, intent(out) :: values(:,:) - character(len=*), intent(in) :: file_name - integer, intent(in) :: n - integer, optional, intent(out) :: rc - - logical :: FileExists - integer :: i,k,face,total, unit, max_procs - integer :: status, N_proc,NF - integer, allocatable :: values_tmp(:), values_(:,:) - - - N_proc = n*6 ! it has been devided by 6. get back the original NY - allocate(values_tmp(N_proc), stat=status) ! no point in checking status - _VERIFY(status) - - inquire(FILE = trim(file_name), EXIST=FileExists) - if ( .not. FileExists) then - print*, file_name // " does not exist" - _RETURN(_FAILURE) - - elseif (MAPL_AM_I_Root(VM)) then - - UNIT = GETFILE ( trim(file_name), form="formatted", rc=status ) - _VERIFY(STATUS) - read(UNIT,*) total, max_procs - if (total /= N_proc) then - print*, "n /= total" - _RETURN(_FAILURE) - endif - do i = 1,total - read(UNIT,*) values_tmp(i) - enddo - call FREE_FILE(UNIT) - endif - - call MAPL_CommsBcast(VM, max_procs, n=1, ROOT=MAPL_Root, rc=status) - call MAPL_CommsBcast(VM, values_tmp, n=N_proc, ROOT=MAPL_Root, rc=status) - _VERIFY(STATUS) - - ! distributed to 6 faces - allocate(values_(max_procs,6)) - values_ = 0 - k = 1 - do NF = 1, 6 - face = 0 - do i = 1, max_procs - values_(i,NF) = values_tmp(k) - face = face + values_tmp(k) - k = k+1 - if (face == this%im_world) exit - enddo - enddo - values = values_ - - _RETURN(_SUCCESS) - - end subroutine get_jms_from_file - - subroutine get_bounds(bounds, label, rc) - type(RealMinMax), intent(out) :: bounds - character(len=*) :: label - integer, optional, intent(out) :: rc - - integer :: i - integer :: n - integer :: status - - call ESMF_ConfigFindLabel(config, label=prefix//label,rc=status) - if (status /= _SUCCESS) then - _RETURN(_SUCCESS) - end if - - ! Must be 2 values: min and max - call ESMF_ConfigGetAttribute(config, bounds%min, rc=status) - _VERIFY(status) - call ESMF_ConfigGetAttribute(config, bounds%max, rc=status) - _VERIFY(status) - - _RETURN(_SUCCESS) - - end subroutine get_bounds - - - end subroutine initialize_from_config_with_prefix - - subroutine halo_init(this, halo_width) - use CubeHaloMod, only: mpp_domain_decomp - class (CubedSphereGridFactory), intent(inout) :: this - integer, optional, intent(in) :: halo_width - integer :: npx, npy - integer :: ng, nregions, grid_type - - nregions = 6 - npx = this%im_world + 1 - npy = npx - ng = 1 ! halowidth - if(present(halo_width)) ng = halo_width - grid_type = 0 ! cubed-sphere - - call mpp_domain_decomp(this%domain, npx, npy, nregions, ng, grid_type, & - & this%nx*this%ny, this%nx, this%ny) - - end subroutine halo_init - - function to_string(this) result(string) - character(len=:), allocatable :: string - class (CubedSphereGridFactory), intent(in) :: this - - _UNUSED_DUMMY(this) - string = 'CubedSphereGridFactory' - - end function to_string - - - - subroutine check_and_fill_consistency(this, unusable, rc) - use MAPL_BaseMod, only: MAPL_DecomposeDim - class (CubedSphereGridFactory), intent(inout) :: this - class (KeywordEnforcer), optional, intent(in) :: unusable - integer, optional, intent(out) :: rc - - integer :: status - character(len=*), parameter :: Iam = MOD_NAME // 'check_and_fill_consistency' - - - if (.not. allocated(this%grid_name)) then - this%grid_name = GRID_NAME_DEFAULT - end if - - if (this%grid_type == UNDEFINED_INTEGER) then - this%grid_type = FV_GRID_TYPE_DEFAULT ! fv default - end if - - ! Check decomposition/bounds - ! WY notes: not necessary for this assert - !_ASSERT(allocated(this%ims) .eqv. allocated(this%jms)) - call verify(this%nx, this%im_world, this%ims, rc=status) - if (allocated(this%jms_2d)) then - _ASSERT(size(this%jms_2d,2)==6) - _ASSERT(sum(this%jms_2d) == 6*this%im_world) - else - call verify(this%ny, this%im_world, this%jms, rc=status) - endif - - _RETURN(_SUCCESS) - - contains - - subroutine verify(n, m_world, ms, rc) - integer, intent(inout) :: n - integer, intent(inout) :: m_world - integer, allocatable, intent(inout) :: ms(:) - integer, optional, intent(out) :: rc - - integer :: status - - if (allocated(ms)) then - _ASSERT(size(ms) > 0) - - if (n == UNDEFINED_INTEGER) then - n = size(ms) - else - _ASSERT(n == size(ms)) - end if - - if (m_world == UNDEFINED_INTEGER) then - m_world = sum(ms) - else - _ASSERT(m_world == sum(ms)) - end if - - else - - _ASSERT(n /= UNDEFINED_INTEGER) - _ASSERT(m_world /= UNDEFINED_INTEGER) - allocate(ms(n), stat=status) - _VERIFY(status) - - call MAPL_DecomposeDim (m_world, ms, n, symmetric=.true.) - - end if - - _RETURN(_SUCCESS) - - end subroutine verify - - end subroutine check_and_fill_consistency - - - elemental subroutine set_with_default_integer(to, from, default) - integer, intent(out) :: to - integer, optional, intent(in) :: from - integer, intent(in) :: default - - if (present(from)) then - to = from - else - to = default - end if - - end subroutine set_with_default_integer - - - elemental subroutine set_with_default_real(to, from, default) - real, intent(out) :: to - real, optional, intent(in) :: from - real, intent(in) :: default - - if (present(from)) then - to = from - else - to = default - end if - - end subroutine set_with_default_real - - subroutine set_with_default_character(to, from, default) - character(len=:), allocatable, intent(out) :: to - character(len=*), optional, intent(in) :: from - character(len=*), intent(in) :: default - - if (present(from)) then - to = from - else - to = default - end if - - end subroutine set_with_default_character - - - elemental subroutine set_with_default_bounds(to, from, default) - type (RealMinMax), intent(out) :: to - type (RealMinMax), optional, intent(in) :: from - type (RealMinMax), intent(in) :: default - - if (present(from)) then - to = from - else - to = default - end if - - end subroutine set_with_default_bounds - - - logical function equals(a, b) - class (CubedSphereGridFactory), intent(in) :: a - class (AbstractGridFactory), intent(in) :: b - - select type (b) - class default - equals = .false. - return - class is (CubedSphereGridFactory) - equals = .true. - - equals = (a%im_world == b%im_world) - if (.not. equals) return - - equals = (a%lm == b%lm) - if (.not. equals) return - - ! same decomposition - equals = all(a%ims == b%ims) - if (.not. equals) return - - if ( allocated(a%jms) .and. allocated(b%jms)) then - equals = all(a%jms == b%jms) - if (.not. equals) return - else - equals = all(a%jms_2d == b%jms_2d) - if (.not. equals) return - endif - - end select - - end function equals - - subroutine initialize_from_esmf_distGrid(this, dist_grid, lon_array, lat_array, unusable, rc) - class (CubedSphereGridFactory), intent(inout) :: this - type (ESMF_DistGrid), intent(in) :: dist_grid - type (ESMF_LocalArray), intent(in) :: lon_array - type (ESMF_LocalArray), intent(in) :: lat_array - class (KeywordEnforcer), optional, intent(in) :: unusable - integer, optional, intent(out) :: rc - - integer :: status - character(len=*), parameter :: Iam = MOD_NAME // 'CubedSphereGridFactory_initialize_from_esmf_distGrid' - - ! not implemented - _ASSERT(.false.) - - end subroutine initialize_from_esmf_distGrid - - subroutine halo(this, array, unusable, halo_width, rc) - use, intrinsic :: iso_fortran_env, only: REAL32 - class (CubedSphereGridFactory), intent(inout) :: this - real(kind=REAL32), intent(inout) :: array(:,:) - class (KeywordEnforcer), optional, intent(in) :: unusable - integer, optional, intent(in) :: halo_width - integer, optional, intent(out) :: rc - - integer :: status - character(len=*), parameter :: Iam = MOD_NAME // 'halo' - - if (.not. this%halo_initialized) then - call this%halo_init(halo_width = halo_width) - this%halo_initialized = .true. - end if - - call mpp_update_domains(array, this%domain) - - _RETURN(_SUCCESS) - - end subroutine halo - - function generate_grid_name(this) result(name) - class (CubedSphereGridFactory), intent(in) :: this - character(len=:), allocatable :: name - - character(len=4) :: im_string - - write(im_string,'(i4.4)') this%im_world - - name = 'CF' // im_string //'x6C' - - end function generate_grid_name - - subroutine a2d2c_factory(this, u, v, lm, getC, unusable, rc) - class (CubedSphereGridFactory), intent(in) :: this - real, intent(inout) :: u(:,:,:) - real, intent(inout) :: v(:,:,:) - integer, intent(in) :: lm - logical, intent(in) :: getC - class (KeywordEnforcer), optional, intent(in) :: unusable - integer, optional, intent(out) :: rc - - character(len=*), parameter :: Iam= MOD_NAME // 'A2D2C' - integer :: status - - ! From GetWeights - call a2d2c(u,v,lm,getC) - - _RETURN(_SUCCESS) - - end subroutine A2D2C_FACTORY - -end module CubedSphereGridFactoryMod diff --git a/DynCore_GridCompMod.F90 b/DynCore_GridCompMod.F90 index a1e27bf..80b31f4 100644 --- a/DynCore_GridCompMod.F90 +++ b/DynCore_GridCompMod.F90 @@ -18,9 +18,8 @@ Module FVdycoreCubed_GridComp ! !USES: use ESMF ! ESMF base class - use MAPL_Mod ! GEOS base class + use MAPL ! GEOS base class use m_set_eta, only: set_eta - use m_chars, only: uppercase ! FV Specific Module use fv_arrays_mod, only: REAL4, REAL8, FVPRC @@ -56,10 +55,7 @@ Module FVdycoreCubed_GridComp fv_getUpdraftHelicity, & ADIABATIC, SW_DYNAMICS, AdvCore_Advection use m_topo_remap, only: dyn_topo_remap - use MAPL_GridManagerMod - use MAPL_RegridderManagerMod - use MAPL_AbstractRegridderMod - use MAPL_RegridderSpecMod + use CubeGridPrototype, only: register_grid_and_regridders ! !PUBLIC MEMBER FUNCTIONS: @@ -2891,6 +2887,7 @@ subroutine Run(gc, import, export, clock, rc) character(len=ESMF_MAXSTR), allocatable :: xlist(:) character(len=ESMF_MAXSTR), allocatable :: biggerlist(:) integer, parameter :: XLIST_MAX = 60 + logical :: isPresent Iam = "Run" call ESMF_GridCompGet( GC, name=COMP_NAME, CONFIG=CF, grid=ESMFGRID, RC=STATUS ) @@ -3097,8 +3094,9 @@ subroutine Run(gc, import, export, clock, rc) firstime=.false. ! get the list of excluded tracers from resource n = 0 - call ESMF_ConfigFindLabel ( CF,'EXCLUDE_ADVECTION_TRACERS_LIST:',rc=STATUS ) - if(STATUS==ESMF_SUCCESS) then + call ESMF_ConfigFindLabel ( CF,'EXCLUDE_ADVECTION_TRACERS_LIST:',isPresent=isPresent,rc=STATUS ) + VERIFY_(STATUS) + if(isPresent) then tend = .false. allocate(xlist(XLIST_MAX), stat=status) @@ -3412,12 +3410,12 @@ subroutine Run(gc, import, export, clock, rc) if ( (qqq%is_r4) .and. associated(qqq%content_r4) ) then if (size(qv)==size(qqq%content_r4)) then qv = qqq%content_r4 - ASSERT_(all(qv >= 0.0)) + _ASSERT(all(qv >= 0.0),'needs informative message') endif elseif (associated(qqq%content)) then if (size(qv)==size(qqq%content)) then qv = qqq%content - ASSERT_(all(qv >= 0.0)) + _ASSERT(all(qv >= 0.0),'needs informative message') endif endif endif @@ -3628,9 +3626,9 @@ subroutine Run(gc, import, export, clock, rc) call MAPL_GetResource(MAPL, ANA_IS_WEIGHTED, Label="ANA_IS_WEIGHTED:", default='NO', RC=STATUS) VERIFY_(STATUS) - ANA_IS_WEIGHTED = uppercase(ANA_IS_WEIGHTED) + ANA_IS_WEIGHTED = ESMF_UtilStringUpperCase(ANA_IS_WEIGHTED) IS_WEIGHTED = adjustl(ANA_IS_WEIGHTED)=="YES" .or. adjustl(ANA_IS_WEIGHTED)=="NO" - ASSERT_( IS_WEIGHTED ) + _ASSERT( IS_WEIGHTED ,'needs informative message') IS_WEIGHTED = adjustl(ANA_IS_WEIGHTED)=="YES" ! Add Analysis Tendencies @@ -6779,7 +6777,7 @@ subroutine ADD_INCS ( STATE,IMPORT,DT,IS_WEIGHTED,RC ) if (TRIM(state%vars%tracer(n)%tname) == 'QILS' ) nwat_tracers = nwat_tracers + 1 enddo ! We must have these first 5 at a minimum - ASSERT_(nwat_tracers == 5) + _ASSERT(nwat_tracers == 5, 'needs informative message') ! Check for QRAIN, QSNOW, QGRAUPEL do n=1,STATE%GRID%NQ if (TRIM(state%vars%tracer(n)%tname) == 'QRAIN' ) nwat_tracers = nwat_tracers + 1 @@ -7338,7 +7336,7 @@ subroutine VertInterp(v2,v3,ple,pp,rc) km = size(ple,3)-1 edge = size(v3,3)==km+1 - ASSERT_(edge .or. size(v3,3)==km) + _ASSERT(edge .or. size(v3,3)==km,'needs informative message') v2 = MAPL_UNDEF @@ -7447,6 +7445,8 @@ subroutine Coldstart(gc, import, export, clock, rc) type (DynState), pointer :: STATE type (DynGrid), pointer :: GRID + logical :: isPresent + ! Tracer Stuff real(r4), pointer :: TRACER(:,:,:) real(REAL8), allocatable :: Q5(:,:,:) @@ -7554,8 +7554,9 @@ subroutine Coldstart(gc, import, export, clock, rc) U(IS:IE,JS:JE,KE) = .001*abs(lats(:,:)) V = 0.0 - call ESMF_ConfigFindLabel( cf, 'AK:', rc = status ) - if (STATUS == 0) then + call ESMF_ConfigFindLabel( cf, 'AK:', isPresent=isPresent, rc = status ) + VERIFY_(STATUS) + if (isPresent) then do L = 0, SIZE(AK)-1 call ESMF_ConfigNextLine ( CF, rc=STATUS ) call ESMF_ConfigGetAttribute( cf, AK(L), rc = status ) @@ -7565,8 +7566,9 @@ subroutine Coldstart(gc, import, export, clock, rc) ak_is_missing = .true. endif - call ESMF_ConfigFindLabel( cf, 'BK:', rc = status ) - if (STATUS == 0) then + call ESMF_ConfigFindLabel( cf, 'BK:', isPresent=isPresent, rc = status ) + VERIFY_(STATUS) + if (isPresent) then do L = 0, SIZE(bk)-1 call ESMF_ConfigNextLine ( CF, rc=STATUS ) call ESMF_ConfigGetAttribute( cf, BK(L), rc = status ) @@ -7578,7 +7580,7 @@ subroutine Coldstart(gc, import, export, clock, rc) if (ak_is_missing .or. bk_is_missing) call set_eta(km, ls, ptop, pint, AK, BK) - ASSERT_(ANY(AK /= 0.0) .or. ANY(BK /= 0.0)) + _ASSERT(ANY(AK /= 0.0) .or. ANY(BK /= 0.0),'needs informative message') do L=lbound(PE,3),ubound(PE,3) PE(:,:,L) = AK(L) + BK(L)*MAPL_P00 enddo @@ -8604,29 +8606,4 @@ function R4_TO_R8(dbl_var) enddo end function - subroutine register_grid_and_regridders() - use MAPL_GridManagerMod, only: grid_manager - use CubedSphereGridFactoryMod, only: CubedSphereGridFactory - use MAPL_RegridderManagerMod, only: regridder_manager - use MAPL_RegridderSpecMod, only: REGRID_METHOD_BILINEAR - use LatLonToCubeRegridderMod - use CubeToLatLonRegridderMod - use CubeToCubeRegridderMod - - type (CubedSphereGridFactory) :: factory - - type (CubeToLatLonRegridder) :: cube_to_latlon_prototype - type (LatLonToCubeRegridder) :: latlon_to_cube_prototype - type (CubeToCubeRegridder) :: cube_to_cube_prototype - - call grid_manager%add_prototype('Cubed-Sphere',factory) - associate (method => REGRID_METHOD_BILINEAR, mgr => regridder_manager) - call mgr%add_prototype('Cubed-Sphere', 'LatLon', method, cube_to_latlon_prototype) - call mgr%add_prototype('LatLon', 'Cubed-Sphere', method, latlon_to_cube_prototype) - call mgr%add_prototype('Cubed-Sphere', 'Cubed-Sphere', method, cube_to_cube_prototype) - end associate - - end subroutine register_grid_and_regridders - - end module FVdycoreCubed_GridComp diff --git a/ESMF_GenerateCSGridDescription.F90 b/ESMF_GenerateCSGridDescription.F90 new file mode 100644 index 0000000..e878291 --- /dev/null +++ b/ESMF_GenerateCSGridDescription.F90 @@ -0,0 +1,949 @@ +#define I_AM_MAIN +#define VERIFY_(A) if(MAPL_VRFY(A,Iam,__LINE__,RC))call MAPL_Abort +!#include "MAPL_Generic.h" + + program ESMF_GenerateCSGridDescription + +! ESMF Framework module + use NetCDF + use ESMF + use MAPL + use fv_grid_utils_mod, only: gnomonic_grids, cell_center2,get_area, direct_transform + use fv_grid_tools_mod, only: mirror_grid + use, intrinsic :: iso_fortran_env, only: REAL64,REAL32 + + implicit none + +! Local variables to be used in the Regrid method calls. +! The code creating and filling these variables is not included in the +! example documentation because those interfaces are not specific to +! Regrid. + type(ESMF_Grid) :: dstgrid + integer :: rc + +!EOC + + real(REAL64), parameter :: PI = 3.14159265358979323846 + character(len=ESMF_MAXSTR) :: IAm + integer :: npets, localPet + integer :: i, j, k + type(ESMF_VM) :: vm + real (ESMF_KIND_R8), dimension(2) :: mincoord + real (ESMF_KIND_R8) :: deltaX, deltaY + integer :: IM_World, JM_World, scrip_size + integer, parameter :: grid_type = 0 + integer, parameter :: KM_WORLD=1 + integer, parameter :: NX=1 + integer, parameter :: NY=6 + integer, parameter :: ntiles=6 + integer, parameter :: ndims=2 + integer :: N + integer :: info + integer :: start(2), cnt(2), hull_num, hull(4) + integer :: UNIT + real(ESMF_KIND_R8), allocatable :: grid_global(:,:,:,:) + real(ESMF_KIND_R8), allocatable :: SCRIP_CenterLat(:), SCRIP_CenterLon(:) + real(ESMF_KIND_R8), allocatable :: SCRIP_CornerLat(:,:), SCRIP_CornerLon(:,:) + real(ESMF_KIND_R8), allocatable :: SCRIP_Area(:) + real(ESMF_KIND_R8) :: alocs(2), node_xy(2,4), node_xy_tmp(2,4), lon_e, lon_w + integer :: cornerdim, gridsize, griddim, rankdim, mask + integer :: cornerlon, cornerlat, centerlon, centerlat,cellarea + integer, allocatable :: IMS(:), JMS(:), sendData(:), GlobalCounts(:), recvCounts(:), recvOffsets(:) + integer, allocatable :: grid_imask(:) + character(len=ESMF_MAXSTR) :: gridname, FMT, FMTIM, FMTJM, FMTPET, title + integer :: NPX, NPY + integer :: isg, ieg + integer :: jsg, jeg + integer :: is, ie + integer :: js, je + integer :: myTile + integer :: npts, tmp, mpiC, mpiC2 + integer :: IG, JG + logical :: do_schmidt + real(ESMF_KIND_R8) :: target_lon, target_lat, stretch_fac + type(ESMF_Config) :: CF + integer :: status + +#include "mpif.h" + + Iam = "ESMF_GenerateCSGridDescription" + + call ESMF_Initialize(logKindFlag=ESMF_LOGKIND_NONE,rc=status) + VERIFY_(status) + + call ESMF_VMGetGlobal(vm, rc=status) + VERIFY_(status) + +! Get number of PETs we are running with +! -------------------------------------- + call ESMF_VMGet(vm, localPet=localPet, petCount=npets, mpiCommunicator=mpiC, rc=status) + VERIFY_(status) + call MPI_Comm_dup(mpic, mpic2, status) + VERIFY_(status) +! Duplicate the MPI communicator not to interfere with ESMF communications. +! The duplicate MPI communicator can be used in any MPI call in the user +! code. Here the MPI_Barrier() routine is called. + call MPI_Barrier(mpic2, status) + VERIFY_(status) + + CF = ESMF_ConfigCreate(rc=status) + VERIFY_(STATUS) + call ESMF_ConfigLoadFile(CF,filename='GenScrip.rc',rc=status) + VERIFY_(STATUS) + call ESMF_ConfigGetAttribute(CF, IM_World, Label='CUBE_DIM:',rc=status) + VERIFY_(STATUS) + JM_WORLD = 6 * IM_WORLD + call ESMF_ConfigGetAttribute(CF, do_schmidt, Label='DO_SCHMIDT:',Default=.false.,rc=status) + VERIFY_(STATUS) + call ESMF_ConfigGetAttribute(CF, stretch_fac, Label='STRETCH_FAC:',rc=status) + if (status /=0) then + if (do_schmidt) then + write(*,*)"Asking for stretch grid without supplying stretch factor" + call MPI_Abort(mpiC,status) + end if + end if + call ESMF_ConfigGetAttribute(CF, target_lon, Label='TARGET_LON:',rc=status) + if (status /=0) then + if (do_schmidt) then + write(*,*)"Asking for stretch grid without supplying target lon" + call MPI_Abort(mpiC,status) + end if + end if + call ESMF_ConfigGetAttribute(CF, target_lat, Label='TARGET_LAT:',rc=status) + if (status /=0) then + if (do_schmidt) then + write(*,*)"Asking for stretch grid without supplying target lat" + call MPI_Abort(mpiC,status) + end if + end if + +! Create destination Cubed-Sphere grid +! ------------------------------------ + allocate( IMS(0:NX-1) ) + allocate( JMS(0:NY-1) ) + + call DecomposeDim_ ( IM_WORLD , IMS , NX ) + call DecomposeDim_ ( JM_WORLD/6, JMS(0:NY/6 -1) , NY/6 ) + do n=2,6 + JMS((n-1)*NY/6 : n*NY/6 -1) = JMS(0:NY/6 -1) + enddo + + deltaX = 2.0*PI/IM_WORLD + deltaY = PI/(JM_WORLD ) + minCoord(1) = 0.0 + minCoord(2) = -PI/2 + + call GET_INT_FORMAT_(IM_WORLD, FMTIM) + call GET_INT_FORMAT_(JM_WORLD, FMTJM) + FMT = '(A,' // trim(FMTIM) //',A,' // trim(FMTJM) // ',A)' + write(gridname,trim(FMT)) 'PE',IM_WORLD,'x',JM_WORLD,'-CF' + dstgrid = ESMF_GridCreate( & + name=gridname, & + countsPerDEDim1=ims, & + countsPerDEDim2=jms, & + indexFlag = ESMF_INDEX_GLOBAL, & + coordDep1 = (/1,2/), & + coordDep2 = (/1,2/), & + rc=status) + VERIFY_(status) + +! Allocate coords at default stagger location +! ------------------------------------------- + call ESMF_GridAddCoord(dstgrid, rc=status) + VERIFY_(status) + + call ESMF_AttributeSet(dstgrid, 'GridType', 'Cubed-Sphere', rc=status) + VERIFY_(status) + + call ESMF_GRID_INTERIOR(dstgrid,isg,ieg,jsg,jeg) +! print*, 'ESMF_GRID_INTERIOR: ', isg, ieg, jsg, jeg + + npx = IM_WORLD + npy = JM_WORLD + myTile = jsg/(npy/ntiles) + + is = isg + ie = ieg + js = jsg - myTile*(npy/ntiles) + je = jeg - myTile*(npy/ntiles) + npts = (npy/ntiles) + if (npts /= npx) then + print*, 'Error npts /= npx', npts, npx + status=1 + endif + VERIFY_(status) + + print*, 'AppGridCreate: ', myTile, is, ie, js, je, npts + + allocate( grid_global(npts+1,npts+1,ndims,ntiles) ) + grid_global=Z'7FFC000000000000' + + call gnomonic_grids(grid_type, npts, grid_global(:,:,1,1), grid_global(:,:,2,1)) + +! mirror_grid assumes that the tile=1 is centered on equator and greenwich meridian Lon[-pi,pi] +! --------------------------------------------------------------------------------------------- + call mirror_grid(grid_global, 0, npts+1, npts+1, 2, 6) + + do n=1,ntiles + do j=1,npts+1 + do i=1,npts+1 + +! Shift the corner away from Japan close to east coast of China +! ------------------------------------------------------------- + if (.not. do_schmidt) grid_global(i,j,1,n) = grid_global(i,j,1,n) - pi/18. + if ( grid_global(i,j,1,n) < 0. ) & + grid_global(i,j,1,n) = grid_global(i,j,1,n) + 2.*pi + if (ABS(grid_global(i,j,1,n)) < 1.e-10) grid_global(i,j,1,n) = 0.0 + if (ABS(grid_global(i,j,2,n)) < 1.e-10) grid_global(i,j,2,n) = 0.0 + enddo + enddo + enddo + +! Clean Up Corners +! ---------------- + grid_global( 1,1:npts+1,:,2)=grid_global(npts+1,1:npts+1,:,1) + grid_global( 1,1:npts+1,:,3)=grid_global(npts+1:1:-1,npts+1,:,1) + grid_global(1:npts+1,npts+1,:,5)=grid_global(1,npts+1:1:-1,:,1) + grid_global(1:npts+1,npts+1,:,6)=grid_global(1:npts+1,1,:,1) + grid_global(1:npts+1, 1,:,3)=grid_global(1:npts+1,npts+1,:,2) + grid_global(1:npts+1, 1,:,4)=grid_global(npts+1,npts+1:1:-1,:,2) + grid_global(npts+1,1:npts+1,:,6)=grid_global(npts+1:1:-1,1,:,2) + grid_global( 1,1:npts+1,:,4)=grid_global(npts+1,1:npts+1,:,3) + grid_global( 1,1:npts+1,:,5)=grid_global(npts+1:1:-1,npts+1,:,3) + grid_global(npts+1,1:npts+1,:,3)=grid_global(1,1:npts+1,:,4) + grid_global(1:npts+1, 1,:,5)=grid_global(1:npts+1,npts+1,:,4) + grid_global(1:npts+1, 1,:,6)=grid_global(npts+1,npts+1:1:-1,:,4) + grid_global( 1,1:npts+1,:,6)=grid_global(npts+1,1:npts+1,:,5) + + !------------------------ + ! Schmidt transformation: + !------------------------ + if ( do_schmidt ) then + target_lon = target_lon*PI/180._8 + target_lat = target_lat*PI/180._8 + do n=1,ntiles + call direct_transform(stretch_fac, 1, npts+1, 1, npts+1, target_lon, target_lat, & + n, grid_global(1:npts+1,1:npts+1,1,n), grid_global(1:npts+1,1:npts+1,2,n)) + enddo + endif + + tmp = (ieg-isg+1)*(jeg-jsg+1) + allocate(SCRIP_CenterLat(tmp),stat=status) + VERIFY_(status) + SCRIP_CenterLat=Z'7FFC000000000000' + allocate(SCRIP_CenterLon(tmp),stat=status) + VERIFY_(status) + SCRIP_CenterLon=Z'7FFC000000000000' + allocate(SCRIP_CornerLat(4,tmp),stat=status) + VERIFY_(status) + SCRIP_CornerLat=Z'7FFC000000000000' + allocate(SCRIP_CornerLon(4,tmp),stat=status) + VERIFY_(status) + SCRIP_CornerLon=Z'7FFC000000000000' + allocate(SCRIP_Area(tmp),stat=status) + VERIFY_(status) + SCRIP_Area=Z'7FFC000000000000' + +! write a separate GMT format multi-segment polygon file from each process +! can be used later to inspect the grid http://gmt.soest.hawaii.edu/ +! ------------------------------------------------------------------------ + call GET_INT_FORMAT_(localPet, FMTPET) + !FMT = '(A,'// trim(FMTIM) //',A,' // trim(FMTPET) //',A)' + !write(filename,trim(FMT)) 'c', IM_WORLD,'.',localPet,'.gmt' + !UNIT=28+localPet + !open(UNIT,file=trim(filename),form='formatted',status='new') + n=1 + do jg=jsg,jeg + do ig=isg,ieg + i=ig + j=jg-myTile*npts + call cell_center2(grid_global(i ,j ,1:2,myTile+1), grid_global(i+1,j ,1:2,myTile+1), & + grid_global(i+1,j+1,1:2,myTile+1), grid_global(i ,j+1,1:2,myTile+1), & + alocs) + SCRIP_CenterLon(n) = alocs(1)*(180._8/PI) + SCRIP_CenterLat(n) = alocs(2)*(180._8/PI) + + node_xy(1,1:4) = (/grid_global(i ,j,1,myTile+1),grid_global(i+1,j,1,myTile+1),grid_global(i+1,j+1,1,myTile+1),grid_global(i,j+1,1,myTile+1)/) + node_xy(2,1:4) = (/grid_global(i ,j,2,myTile+1),grid_global(i+1,j,2,myTile+1),grid_global(i+1,j+1,2,myTile+1),grid_global(i,j+1,2,myTile+1)/) + node_xy_tmp = node_xy + +! Correct for the periodic boundary at 0/360 +! ------------------------------------------ + lon_w = min( grid_global(i ,j,1,myTile+1),grid_global(i+1,j,1,myTile+1),grid_global(i+1,j+1,1,myTile+1),grid_global(i,j+1,1,myTile+1) ) + lon_e = max( grid_global(i ,j,1,myTile+1),grid_global(i+1,j,1,myTile+1),grid_global(i+1,j+1,1,myTile+1),grid_global(i,j+1,1,myTile+1) ) + if ( abs(lon_e - lon_w) > 1.5_8*pi .and. (SCRIP_CenterLon(n) < pi) ) then + where(node_xy(1,:) > pi) node_xy_tmp(1,:) = node_xy(1,:) - 2._8*pi + elseif ( abs(lon_e - lon_w) > 1.5_8*pi .and. (SCRIP_CenterLon(n) > pi) ) then + where(node_xy(1,:) < pi) node_xy_tmp(1,:) = node_xy(1,:) + 2._8*pi + endif + call points_hull_2d(4, node_xy_tmp, hull_num, hull) + if(ANY(hull==0)) then + write(*,100)'Zero Hull ', grid_global(i ,j,1,myTile+1),grid_global(i+1,j,1,myTile+1),grid_global(i+1,j+1,1,myTile+1),grid_global(i,j+1,1,myTile+1) + write(*,100)'Zero Hull ', node_xy_tmp(1,:) + endif + !write(UNIT,103) '>' + do k=1,4 + SCRIP_CornerLon(k,n) = node_xy(1,hull(k))*(180._8/PI) + SCRIP_CornerLat(k,n) = node_xy(2,hull(k))*(180._8/PI) + !write(UNIT,102) SCRIP_CornerLon(k,n), SCRIP_CornerLat(k,n) + enddo + !SCRIP_Area(n) = get_area(node_xy(:,hull(1)),node_xy(:,hull(2)),node_xy(:,hull(3)),node_xy(:,hull(4)),1.0d0) + SCRIP_Area(n) = get_area(grid_global(i ,j ,1:2,myTile+1), grid_global(i,j+1 ,1:2,myTile+1), & + grid_global(i+1,j,1:2,myTile+1), grid_global(i+1,j+1,1:2,myTile+1),1.0d0) + n=n+1 + enddo + enddo + + 100 format(a,4f20.15) + 101 format(a,f20.15) + 102 format(2f20.15) + 103 format(a) + + deallocate( grid_global ) + deallocate( IMS ) + deallocate( JMS ) + + scrip_size = IM_World*JM_World + call MPI_Info_create(info, status) + VERIFY_(status) + call MPI_Info_set(info, "cb_buffer_size", "1048576", status) + VERIFY_(status) + + status = nf90_create(trim(gridname)//'.nc4', IOR(NF90_MPIIO,IOR(NF90_CLOBBER,NF90_NETCDF4)),unit,comm=mpic2,info=info) + if(status /= nf90_noerr) then + print*,'Error creating file ',status + print*, NF90_STRERROR(status) + stop + endif + + FMT = '(A,' // ',A,' //',A)' + write(title,trim(FMT)) 'GMAO ',trim(gridname),' Grid' + status = nf90_put_att(unit,NF90_GLOBAL, 'title', trim(title)) + if(status /= nf90_noerr) then + print*,'Error setting title',status + print*, NF90_STRERROR(status) + stop + endif + status = nf90_put_att(unit,NF90_GLOBAL, 'GridDescriptionFormat','SCRIP') + if(status /= nf90_noerr) then + print*,'Error setting GridDescriptionFormat',status + print*, NF90_STRERROR(status) + stop + endif + if (do_Schmidt) then + status = nf90_put_att(unit,NF90_GLOBAL, 'TargetLon', target_lon*180._8/PI) + if(status /= nf90_noerr) then + print*,'Error setting TargetLon',status + print*, NF90_STRERROR(status) + stop + endif + status = nf90_put_att(unit,NF90_GLOBAL, 'TargetLat', target_lat*180._8/PI) + if(status /= nf90_noerr) then + print*,'Error setting TargetLat',status + print*, NF90_STRERROR(status) + stop + endif + status = nf90_put_att(unit,NF90_GLOBAL, 'StretchFactor', stretch_fac) + if(status /= nf90_noerr) then + print*,'Error setting StretchFactor',status + print*, NF90_STRERROR(status) + stop + endif + end if + + status = NF90_def_dim(UNIT, 'grid_size', SCRIP_size, gridsize) + if(status /= nf90_noerr) then + print*,'Error defining grid_size',status + print*, NF90_STRERROR(status) + stop + endif + status = nf90_def_dim(unit,'grid_corners', 4, cornerdim) + if(status /= nf90_noerr) then + print*,'Error defining grid_corners',status + print*, NF90_STRERROR(status) + stop + endif + +! Peggy Li suggested setting grid_rank=1 and grid_dims=1 +! so that ESMF will treat the grid as unstructured +! ------------------------------------------------------ + status = NF90_DEF_DIM(UNIT, 'grid_rank' , 1, rankdim) + if(status /= nf90_noerr) then + print*,'Error defining grid_rank',status + print*, NF90_STRERROR(status) + stop + endif + +! Grid dimensions +! --------------- + status = nf90_def_var(UNIT, "grid_dims", NF90_INT, [rankdim], griddim) + if(status /= nf90_noerr) then + print*,'Error defining grid_dims',status + print*, NF90_STRERROR(status) + stop + endif + +! Grid mask +! --------- + status = nf90_def_var(UNIT, "grid_imask", NF90_INT, [gridsize], mask) + if(status /= nf90_noerr) then + print*,'Error defining grid_imask',status + print*, NF90_STRERROR(status) + stop + endif + status = nf90_put_att(unit,mask,"units","unitless") + +! cell center Longitude variable +! ------------------------------ + status = nf90_def_var(UNIT, "grid_center_lon", NF90_DOUBLE, [gridsize], centerlon) + if(status /= nf90_noerr) then + print*,'Error defining cell center lons',status + print*, NF90_STRERROR(status) + stop + endif + status = nf90_put_att(UNIT, centerlon, "units" , "degrees") + +! cell center Latitude variable +! ----------------------------- + status = nf90_def_var(UNIT, "grid_center_lat", NF90_DOUBLE, [gridsize], centerlat) + if(status /= nf90_noerr) then + print*,'Error defining cell center lats',status + print*, NF90_STRERROR(status) + stop + endif + status = nf90_put_att(UNIT, centerlat, "units" , "degrees") + +! cell corner Longitude variable +! ------------------------------ + status = nf90_def_var(UNIT, "grid_corner_lon", NF90_DOUBLE, [cornerdim,gridsize], cornerlon) + if(status /= nf90_noerr) then + print*,'Error defining cell corner lons',status + print*, NF90_STRERROR(status) + stop + endif + status = nf90_put_att(UNIT, cornerlon, "units" , "degrees") + +! cell corner Latitude variable +! ----------------------------- + status = nf90_def_var(UNIT, "grid_corner_lat", NF90_DOUBLE, [cornerdim,gridsize], cornerlat) + if(status /= nf90_noerr) then + print*,'Error defining cell corner lats',status + print*, NF90_STRERROR(status) + stop + endif + status = nf90_put_att(UNIT, cornerlat, "units" ,"degrees") + + status = nf90_def_var(UNIT, "grid_area", NF90_DOUBLE, [gridsize], cellarea) + if(status /= nf90_noerr) then + print*,'Error defining cell area',status + print*, NF90_STRERROR(status) + stop + endif + status = nf90_put_att(UNIT, cellarea, "units" , "radians^2") + + status = nf90_enddef(UNIT) + if(status /= nf90_noerr) then + print*,'Error exitting define mode',status + print*, NF90_STRERROR(status) + stop + endif + + rc = NF90_PUT_VAR(UNIT, griddim, 1) + if(rc /= nf90_noerr) then + print*,'Error writing griddim',rc + print*, NF90_STRERROR(rc) + stop + endif + + allocate (sendData(1),GlobalCounts(npets), recvCounts(npets), recvOffsets(npets), stat=status) + VERIFY_(status) + sendData = tmp + recvCounts=1 + recvOffsets=0 + do i=2, npets + recvOffsets(i) = recvOffsets(i-1) + recvCounts(i-1) + enddo + call ESMF_VMGatherV(vm,sendData=sendData,sendCount=1,recvData=GlobalCounts,recvCounts=recvCounts,recvOffsets=recvOffsets,rootPet=0,rc=status) + VERIFY_(status) + call ESMF_VMBroadcast(vm,bcstData=GlobalCounts,count=npets,rootPet=0, rc=status) + VERIFY_(status) + if (localPet == 0) print*,GlobalCounts + + start=1 + do i=1,localPet + start(1) = start(1)+GlobalCounts(i) + enddo + + cnt(1) = tmp; cnt(2)=1 + status = NF90_PUT_VAR(UNIT, centerlon, SCRIP_CenterLon, start=start, count=cnt) + if(status /= nf90_noerr) then + print*,'Error writing CenterLons ',status + print*, NF90_STRERROR(status) + stop + endif + + status = NF90_PUT_VAR(UNIT, centerlat, SCRIP_CenterLat, start=start, count=cnt) + if(status /= nf90_noerr) then + print*,'Error writing CenterLats ',status + print*, NF90_STRERROR(status) + stop + endif + + status = NF90_PUT_VAR(UNIT, cellarea, SCRIP_Area, start=start, count=cnt) + if(status /= nf90_noerr) then + print*,'Error writing CenterLats ',status + print*, NF90_STRERROR(status) + stop + endif +! Grid mask +! --------- + allocate(grid_imask(tmp), stat=status) + VERIFY_(status) + grid_imask = 1 + status = NF90_PUT_VAR(UNIT, mask, grid_imask, start=start, count=cnt) + if(status /= nf90_noerr) then + print*,'Error writing grid_imask',status + print*, NF90_STRERROR(status) + stop + endif + deallocate(grid_imask) + + start(2)=start(1) + start(1)=1 + cnt(1)=4 + cnt(2)=tmp + + status = NF90_PUT_VAR(UNIT, cornerlat, SCRIP_CornerLat, start=start, count=cnt) + if(status /= nf90_noerr) then + print*,'Error writing CornerLats ',status + print*, NF90_STRERROR(status) + stop + endif + + status = NF90_PUT_VAR(UNIT, cornerlon, SCRIP_CornerLon, start=start, count=cnt) + if(status /= nf90_noerr) then + print*,'Error writing CornerLons ',status + print*, NF90_STRERROR(status) + stop + endif + + status = NF90_CLOSE(UNIT) + if(status /= nf90_noerr) then + print*,'Error closing output file ', status + print*, NF90_STRERROR(status) + stop + endif + call ESMF_VMBarrier(vm, rc=status) + + deallocate(SCRIP_CenterLat) + deallocate(SCRIP_CenterLon) + deallocate(SCRIP_CornerLat) + deallocate(SCRIP_CornerLon) + deallocate(sendData) + deallocate(GlobalCounts) + deallocate(recvCounts) + deallocate(recvOffsets) + + call ESMF_Finalize(rc=status) + + contains + + subroutine GET_INT_FORMAT_(N, FMT) + integer :: N + character(len=*) :: FMT + + IF(N < 10) THEN + FMT = 'I1' + ELSE IF (N< 100) THEN + FMT = 'I2' + ELSE IF (N< 1000) THEN + FMT = 'I3' + ELSE IF (N< 10000) THEN + FMT = 'I4' + else + FMT = 'I5' + end IF + end subroutine GET_INT_FORMAT_ + + subroutine DecomposeDim_( dim_world,dim,NDEs ) +! +! From FMS/MPP +! + implicit none + integer dim_world, NDEs + integer dim(0:NDEs-1) + + integer :: is,ie,isg,ieg + integer :: ndiv,ndivs,imax,ndmax,ndmirror,n + integer :: ibegin(0:NDEs-1) + integer :: iend(0:NDEs-1) + + logical :: symmetrize + logical :: even, odd + even(n) = (mod(n,2).EQ.0) + odd (n) = (mod(n,2).EQ.1) + + isg = 1 + ieg = dim_world + ndivs = NDEs + + is = isg + n = 0 + do ndiv=0,ndivs-1 + !modified for mirror-symmetry + !original line + ! ie = is + CEILING( float(ieg-is+1)/(ndivs-ndiv) ) - 1 + + !problem of dividing nx points into n domains maintaining symmetry + !i.e nx=18 n=4 4554 and 5445 are solutions but 4455 is not. + !this will always work for nx even n even or odd + !this will always work for nx odd, n odd + !this will never work for nx odd, n even: for this case we supersede the mirror calculation + ! symmetrize = .NOT. ( mod(ndivs,2).EQ.0 .AND. mod(ieg-isg+1,2).EQ.1 ) + !nx even n odd fails if n>nx/2 + symmetrize = ( even(ndivs) .AND. even(ieg-isg+1) ) .OR. & + ( odd(ndivs) .AND. odd(ieg-isg+1) ) .OR. & + ( odd(ndivs) .AND. even(ieg-isg+1) .AND. ndivs.LT.(ieg-isg+1)/2 ) + + !mirror domains are stored in the list and retrieved if required. + if( ndiv.EQ.0 )then + !initialize max points and max domains + imax = ieg + ndmax = ndivs + end if + !do bottom half of decomposition, going over the midpoint for odd ndivs + if( ndiv.LT.(ndivs-1)/2+1 )then + !domain is sized by dividing remaining points by remaining domains + ie = is + CEILING( REAL(imax-is+1)/(ndmax-ndiv) ) - 1 + ndmirror = (ndivs-1) - ndiv !mirror domain + if( ndmirror.GT.ndiv .AND. symmetrize )then !only for domains over the midpoint + !mirror extents, the max(,) is to eliminate overlaps + ibegin(ndmirror) = max( isg+ieg-ie, ie+1 ) + iend(ndmirror) = max( isg+ieg-is, ie+1 ) + imax = ibegin(ndmirror) - 1 + ndmax = ndmax - 1 + end if + else + if( symmetrize )then + !do top half of decomposition by retrieving saved values + is = ibegin(ndiv) + ie = iend(ndiv) + else + ie = is + CEILING( REAL(imax-is+1)/(ndmax-ndiv) ) - 1 + end if + end if + dim(ndiv) = ie-is+1 + is = ie + 1 + end do + + end subroutine DecomposeDim_ + +!logical function MAPL_VRFY(A,iam,line,rc) + !integer, intent(IN ) :: A + !character*(*), intent(IN ) :: iam + !integer, intent(IN ) :: line + !integer, optional, intent(OUT) :: RC + !MAPL_VRFY = A/=ESMF_SUCCESS + !if(MAPL_VRFY)then + !if(present(RC)) then + !print'(A40,I10)',Iam,line + !RC=A + !endif + !endif +!end function MAPL_VRFY + + !subroutine MAPL_Abort + !call ESMF_Finalize(endflag = ESMF_END_ABORT) + !end subroutine MAPL_Abort + + !subroutine ESMF_GRID_INTERIOR(GRID,I1,IN,J1,JN) + !type (ESMF_Grid), intent(IN) :: grid + !integer, intent(OUT) :: I1, IN, J1, JN + +!! local vars + !integer :: status + !character(len=ESMF_MAXSTR) :: IAm='ESMF_GridInterior' + + !type (ESMF_DistGrid) :: distGrid + !type(ESMF_DELayout) :: LAYOUT + !integer, allocatable :: AL(:,:) + !integer, allocatable :: AU(:,:) + !integer :: nDEs + !integer :: deId + !integer :: gridRank + !integer :: deList(1) + + !call ESMF_GridGet (GRID, dimCount=gridRank, distGrid=distGrid, rc=STATUS) + !call ESMF_DistGridGet(distGRID, delayout=layout, rc=STATUS) + !call ESMF_DELayoutGet(layout, deCount =nDEs, localDeList=deList, rc=status) + !deId = deList(1) + + !allocate (AL(gridRank,0:nDEs-1), stat=status) + !allocate (AU(gridRank,0:nDEs-1), stat=status) + + !call ESMF_DistGridGet(distgrid, & + !minIndexPDe=AL, maxIndexPDe=AU, rc=status) + + !I1 = AL(1, deId) + !IN = AU(1, deId) +!! ASSERT_(gridRank > 1) !ALT: tilegrid is 1d (without RC this only for info) + !J1 = AL(2, deId) + !JN = AU(2, deId) + !deallocate(AU, AL) + + !end subroutine ESMF_GRID_INTERIOR + +!function angle_rad_2d ( p1, p2, p3 ) +subroutine angle_rad_2d ( p1, p2, p3, res ) +!*****************************************************************************80 +! +!! ANGLE_RAD_2D returns the angle swept out between two rays in 2D. +! +! Discussion: +! +! Except for the zero angle case, it should be true that +! +! ANGLE_RAD_2D ( P1, P2, P3 ) + ANGLE_RAD_2D ( P3, P2, P1 ) = 2 * PI +! +! P1 +! / +! / +! / +! / +! P2--------->P3 +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 15 January 2005 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, real ( REAL64 ) P1(2), P2(2), P3(2), define the rays +! P1 - P2 and P3 - P2 which define the angle. +! +! Output, real ( REAL64 ) ANGLE_RAD_2D, the angle swept out by the rays, +! in radians. 0 <= ANGLE_RAD_2D < 2 * PI. If either ray has zero +! length, then ANGLE_RAD_2D is set to 0. + implicit none + + integer, parameter :: dim_num = 2 + + real (REAL64), parameter :: pi = 3.141592653589793D+00 + real (REAL64) p(dim_num) + real (REAL64) p1(dim_num) + real (REAL64) p2(dim_num) + real (REAL64) p3(dim_num) + real (REAL64) res + + p(1) = ( p3(1) - p2(1) ) * ( p1(1) - p2(1) ) & + + ( p3(2) - p2(2) ) * ( p1(2) - p2(2) ) + + + p(2) = ( p3(1) - p2(1) ) * ( p1(2) - p2(2) ) & + - ( p3(2) - p2(2) ) * ( p1(1) - p2(1) ) + + if ( p(1) == 0.0D+00 .and. p(2) == 0.0D+00 ) then + res = 0.0D+00 + return + end if + + res = atan2 ( p(2), p(1) ) + + if ( res < 0.0D+00 ) then + res = res + 2.0D+00 * pi + end if + + return +end + +subroutine points_hull_2d ( node_num, node_xy, hull_num, hull ) + +!*****************************************************************************80 +! +!! POINTS_HULL_2D computes the convex hull of 2D points. +! +! Discussion: +! +! The work involved is N*log(H), where N is the number of points, and H is +! the number of points that are on the hull. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 12 June 2006 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, integer NODE_NUM, the number of nodes. +! +! Input, real ( REAL64 ) NODE_XY(2,NODE_NUM), the coordinates of the nodes. +! +! Output, integer HULL_NUM, the number of nodes that lie on +! the convex hull. +! +! Output, integer HULL(NODE_NUM). Entries 1 through HULL_NUM +! contain the indices of the nodes that form the convex hull, in order. +! + implicit none + + integer node_num + real (REAL64) angle + real (REAL64) angle_max + real (REAL64) di + real (REAL64) dr + integer first + integer hull(node_num) + integer hull_num + integer i + real (REAL64) node_xy(2,node_num) + real (REAL64) p_xy(2) + integer q + real (REAL64) q_xy(2) + integer r + real (REAL64) r_xy(2) + + if ( node_num < 1 ) then + hull_num = 0 + return + end if +! +! If NODE_NUM = 1, the hull is the point. +! + if ( node_num == 1 ) then + hull_num = 1 + hull(1) = 1 + return + end if +! +! If NODE_NUM = 2, then the convex hull is either the two distinct points, +! or possibly a single (repeated) point. +! + if ( node_num == 2 ) then + + if ( node_xy(1,1) /= node_xy(1,2) .or. node_xy(2,1) /= node_xy(2,2) ) then + hull_num = 2 + hull(1) = 1 + hull(2) = 2 + else + hull_num = 1 + hull(1) = 1 + end if + + return + + end if +! +! Find the leftmost point and call it "Q". +! In case of ties, take the bottom-most. +! + q = 1 + do i = 2, node_num + if ( node_xy(1,i) < node_xy(1,q) .or. & + ( node_xy(1,i) == node_xy(1,q) .and. node_xy(2,i) < node_xy(2,q) ) ) then + q = i + end if + end do + + q_xy(1:2) = node_xy(1:2,q) +! +! Remember the starting point, so we know when to stop! +! + first = q + hull_num = 1 + hull(1) = q +! +! For the first point, make a dummy previous point, 1 unit south, +! and call it "P". +! + p_xy(1) = q_xy(1) + p_xy(2) = q_xy(2) - 1.0D+00 +! +! Now, having old point P, and current point Q, find the new point R +! so the angle PQR is maximal. +! +! Watch out for the possibility that the two nodes are identical. +! + do + + r = 0 + angle_max = 0.0D+00 + + do i = 1, node_num + + if ( i /= q .and. & + ( node_xy(1,i) /= q_xy(1) .or. node_xy(2,i) /= q_xy(2) ) ) then + + call angle_rad_2d(p_xy, q_xy, node_xy(1:2,i),angle) + !angle = angle_rad_2d ( p_xy, q_xy, node_xy(1:2,i) ) + + if ( r == 0 .or. angle_max < angle ) then + + r = i + r_xy(1:2) = node_xy(1:2,r) + angle_max = angle +! +! In case of ties, choose the nearer point. +! + else if ( r /= 0 .and. angle == angle_max ) then + + di = ( node_xy(1,i) - q_xy(1) )**2 + ( node_xy(2,i) - q_xy(2) )**2 + dr = ( r_xy(1) - q_xy(1) )**2 + ( r_xy(2) - q_xy(2) )**2 + + if ( di < dr ) then + r = i + r_xy(1:2) = node_xy(1:2,r) + angle_max = angle + end if + + end if + + end if + + end do +! +! We are done when we have returned to the first point on the convex hull. +! + if ( r == first ) then + exit + end if + + hull_num = hull_num + 1 + if ( node_num < hull_num ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'POINTS_HULL_2D - Fatal error!' + write ( *, '(a)' ) ' The algorithm has failed.' + stop + end if +! +! Add point R to convex hull. +! + hull(hull_num) = r +! +! Set P := Q, Q := R, and prepare to search for next point R. +! + q = r + + p_xy(1:2) = q_xy(1:2) + q_xy(1:2) = r_xy(1:2) + + end do + + return +end + + end program ESMF_GenerateCSGridDescription diff --git a/Externals.cfg b/Externals.cfg index c58bdaa..3a0ace2 100644 --- a/Externals.cfg +++ b/Externals.cfg @@ -2,7 +2,7 @@ required = True repo_url = git@github.com:GEOS-ESM/GFDL_atmos_cubed_sphere.git local_path = ./@fvdycore -tag = geos/v1.0.2 +tag = geos/v1.1.0 protocol = git [externals_description] diff --git a/FV_StateMod.F90 b/FV_StateMod.F90 index fbf07dd..4ca7f32 100644 --- a/FV_StateMod.F90 +++ b/FV_StateMod.F90 @@ -8,7 +8,7 @@ module FV_StateMod ! !USES: #if defined( MAPL_MODE ) use ESMF ! ESMF base class - use MAPL_Mod ! MAPL base class + use MAPL ! MAPL base class #endif use MAPL_ConstantsMod, only: MAPL_CP, MAPL_RGAS, MAPL_RVAP, MAPL_GRAV, MAPL_RADIUS, & @@ -165,8 +165,8 @@ module FV_StateMod ! lat-lon = 1 ! cubed-sphere = 6 - real(REAL8), allocatable :: DXC(:,:) ! local C-Gird DeltaX - real(REAL8), allocatable :: DYC(:,:) ! local C-Gird DeltaY + real(REAL8), allocatable :: DXC(:,:) ! local C-Grid DeltaX + real(REAL8), allocatable :: DYC(:,:) ! local C-Grid DeltaY real(REAL8), allocatable :: AREA(:,:) ! local cell area real(REAL8) :: GLOBALAREA ! global area @@ -579,7 +579,7 @@ subroutine FV_Setup(GC,LAYOUT_FILE, RC) !! Setup GFDL microphysics module call gfdl_cloud_microphys_init() - ASSERT_(DT > 0.0) + _ASSERT(DT > 0.0, 'needs informative message') call WRITE_PARALLEL("Dynamics PE Layout ") call WRITE_PARALLEL(FV_Atm(1)%layout(1) ,format='("NPES_X : ",( I3))') @@ -592,7 +592,7 @@ subroutine FV_Setup(GC,LAYOUT_FILE, RC) ! if needed, we could compute, ks by count(BK==0.0) ! then FV will try to run slightly more efficient code ! So far, GEOS-5 has used ks = 0 - ASSERT_(ks <= FV_Atm(1)%flagstruct%NPZ+1) + _ASSERT(ks <= FV_Atm(1)%flagstruct%NPZ+1,'needs informative message') call WRITE_PARALLEL(ks , & format='("Number of true pressure levels =", I5)' ) @@ -1159,7 +1159,7 @@ subroutine FV_Run (STATE, CLOCK, GC, RC) if (TRIM(state%vars%tracer(n)%tname) == 'QILS' ) nwat_tracers = nwat_tracers + 1 enddo ! We must have these first 5 at a minimum - ASSERT_(nwat_tracers == 5) + _ASSERT(nwat_tracers == 5, 'needs informative message') ! Check for CLLS, CLCN, QRAIN, QSNOW, QGRAUPEL do n=1,STATE%GRID%NQ if (TRIM(state%vars%tracer(n)%tname) == 'CLLS' ) nwat_tracers = nwat_tracers + 1 @@ -1176,7 +1176,7 @@ subroutine FV_Run (STATE, CLOCK, GC, RC) endif endif if (FV_Atm(1)%flagstruct%do_sat_adj) then - ASSERT_(FV_Atm(1)%flagstruct%nwat == 6) + _ASSERT(FV_Atm(1)%flagstruct%nwat == 6, 'needs informative message') endif STATE%VARS%nwat = FV_Atm(1)%flagstruct%nwat endif @@ -1190,7 +1190,7 @@ subroutine FV_Run (STATE, CLOCK, GC, RC) (FV_Atm(1)%flagstruct%nwat == 1) .OR. & (FV_Atm(1)%flagstruct%nwat == 3) .OR. & (FV_Atm(1)%flagstruct%nwat == 6) ) - ASSERT_( NWAT_TEST ) + _ASSERT( NWAT_TEST , 'needs informative message') select case ( FV_Atm(1)%flagstruct%nwat ) case (6) FV_Atm(1)%ncnst = STATE%GRID%NQ + 3 ! NQ + Combined QLIQ,QICE,QCLD @@ -1245,16 +1245,16 @@ subroutine FV_Run (STATE, CLOCK, GC, RC) if (.not. ADIABATIC) then select case ( FV_Atm(1)%flagstruct%nwat ) case (0) - ASSERT_(FV_Atm(1)%ncnst == STATE%GRID%NQ) + _ASSERT(FV_Atm(1)%ncnst == STATE%GRID%NQ, 'needs informative message') case (1) - ASSERT_(FV_Atm(1)%ncnst >= 5) - ASSERT_(FV_Atm(1)%ncnst == STATE%GRID%NQ) + _ASSERT(FV_Atm(1)%ncnst >= 5, 'needs informative message') + _ASSERT(FV_Atm(1)%ncnst == STATE%GRID%NQ, 'needs informative message') case (3) - ASSERT_(FV_Atm(1)%ncnst >= 7) - ASSERT_(FV_Atm(1)%ncnst == STATE%GRID%NQ + 2) + _ASSERT(FV_Atm(1)%ncnst >= 7, 'needs informative message') + _ASSERT(FV_Atm(1)%ncnst == STATE%GRID%NQ + 2, 'needs informative message') case (6) - ASSERT_(FV_Atm(1)%ncnst >= 13) - ASSERT_(FV_Atm(1)%ncnst == STATE%GRID%NQ + 3) + _ASSERT(FV_Atm(1)%ncnst >= 13, 'needs informative message') + _ASSERT(FV_Atm(1)%ncnst == STATE%GRID%NQ + 3, 'needs informative message') end select FV_Atm(1)%q(:,:,:,:) = 0.0 if (FV_Atm(1)%flagstruct%nwat > 0) then @@ -1410,36 +1410,36 @@ subroutine FV_Run (STATE, CLOCK, GC, RC) ! Verify select case (FV_Atm(1)%flagstruct%nwat) case (6) - ASSERT_(nn == 13) ! Q, QLCN, QLLS, QICN, QILS, CLLS, CLCN, QRAIN, QSNOW, QGRAUPEL, QLIQ, QICE, QCLD - ASSERT_(SPHU_FILLED) - ASSERT_(QLIQ_FILLED) - ASSERT_(QICE_FILLED) - ASSERT_(RAIN_FILLED) - ASSERT_(SNOW_FILLED) - ASSERT_(GRPL_FILLED) - ASSERT_(QCLD_FILLED) - ASSERT_(QLCN_FILLED) - ASSERT_(QLLS_FILLED) - ASSERT_(QICN_FILLED) - ASSERT_(QILS_FILLED) - ASSERT_(CLCN_FILLED) - ASSERT_(CLLS_FILLED) + _ASSERT(nn == 13, 'needs informative message') ! Q, QLCN, QLLS, QICN, QILS, CLLS, CLCN, QRAIN, QSNOW, QGRAUPEL, QLIQ, QICE, QCLD + _ASSERT(SPHU_FILLED, 'needs informative message') + _ASSERT(QLIQ_FILLED, 'needs informative message') + _ASSERT(QICE_FILLED, 'needs informative message') + _ASSERT(RAIN_FILLED, 'needs informative message') + _ASSERT(SNOW_FILLED, 'needs informative message') + _ASSERT(GRPL_FILLED, 'needs informative message') + _ASSERT(QCLD_FILLED, 'needs informative message') + _ASSERT(QLCN_FILLED, 'needs informative message') + _ASSERT(QLLS_FILLED, 'needs informative message') + _ASSERT(QICN_FILLED, 'needs informative message') + _ASSERT(QILS_FILLED, 'needs informative message') + _ASSERT(CLCN_FILLED, 'needs informative message') + _ASSERT(CLLS_FILLED, 'needs informative message') case (3) - ASSERT_(nn == 7) ! Q, QLCN, QLLS, QICN, QILS, QLIQ, QICE - ASSERT_(SPHU_FILLED) - ASSERT_(QLIQ_FILLED) - ASSERT_(QICE_FILLED) - ASSERT_(QLCN_FILLED) - ASSERT_(QLLS_FILLED) - ASSERT_(QICN_FILLED) - ASSERT_(QILS_FILLED) + _ASSERT(nn == 7, 'needs informative message') ! Q, QLCN, QLLS, QICN, QILS, QLIQ, QICE + _ASSERT(SPHU_FILLED, 'needs informative message') + _ASSERT(QLIQ_FILLED, 'needs informative message') + _ASSERT(QICE_FILLED, 'needs informative message') + _ASSERT(QLCN_FILLED, 'needs informative message') + _ASSERT(QLLS_FILLED, 'needs informative message') + _ASSERT(QICN_FILLED, 'needs informative message') + _ASSERT(QILS_FILLED, 'needs informative message') case (1) - ASSERT_(nn == 5) ! Q, QLCN, QLLS, QICN, QILS - ASSERT_(SPHU_FILLED) - ASSERT_(QLCN_FILLED) - ASSERT_(QLLS_FILLED) - ASSERT_(QICN_FILLED) - ASSERT_(QILS_FILLED) + _ASSERT(nn == 5, 'needs informative message') ! Q, QLCN, QLLS, QICN, QILS + _ASSERT(SPHU_FILLED, 'needs informative message') + _ASSERT(QLCN_FILLED, 'needs informative message') + _ASSERT(QLLS_FILLED, 'needs informative message') + _ASSERT(QICN_FILLED, 'needs informative message') + _ASSERT(QILS_FILLED, 'needs informative message') end select endif !nwat > 0 select case (FV_Atm(1)%flagstruct%nwat) @@ -1503,7 +1503,7 @@ subroutine FV_Run (STATE, CLOCK, GC, RC) endif enddo end select - ASSERT_(nn == FV_Atm(1)%ncnst) + _ASSERT(nn == FV_Atm(1)%ncnst, 'needs informative message') else if (mpp_pe()==0) print*, 'Running In Adiabatic Mode' do n=1,STATE%GRID%NQ @@ -1514,7 +1514,7 @@ subroutine FV_Run (STATE, CLOCK, GC, RC) FV_Atm(1)%q(isc:iec,jsc:jec,1:npz,nn) = state%vars%tracer(n)%content(:,:,:) endif enddo - ASSERT_(nn == FV_Atm(1)%ncnst) + _ASSERT(nn == FV_Atm(1)%ncnst, 'needs informative message') endif myDT = state%dt @@ -1866,36 +1866,36 @@ subroutine FV_Run (STATE, CLOCK, GC, RC) ! Verify select case (FV_Atm(1)%flagstruct%nwat) case (6) - ASSERT_(nn == 13) ! Q, QLCN, QLLS, QICN, QILS, CLLS, CLCN, QRAIN, QSNOW, QGRAUPEL, QLIQ, QICE, QCLD - ASSERT_(SPHU_FILLED) - ASSERT_(QLIQ_FILLED) - ASSERT_(QICE_FILLED) - ASSERT_(RAIN_FILLED) - ASSERT_(SNOW_FILLED) - ASSERT_(GRPL_FILLED) - ASSERT_(QCLD_FILLED) - ASSERT_(QLCN_FILLED) - ASSERT_(QLLS_FILLED) - ASSERT_(QICN_FILLED) - ASSERT_(QILS_FILLED) - ASSERT_(CLCN_FILLED) - ASSERT_(CLLS_FILLED) + _ASSERT(nn == 13, 'needs informative message') ! Q, QLCN, QLLS, QICN, QILS, CLLS, CLCN, QRAIN, QSNOW, QGRAUPEL, QLIQ, QICE, QCLD + _ASSERT(SPHU_FILLED, 'needs informative message') + _ASSERT(QLIQ_FILLED, 'needs informative message') + _ASSERT(QICE_FILLED, 'needs informative message') + _ASSERT(RAIN_FILLED, 'needs informative message') + _ASSERT(SNOW_FILLED, 'needs informative message') + _ASSERT(GRPL_FILLED, 'needs informative message') + _ASSERT(QCLD_FILLED, 'needs informative message') + _ASSERT(QLCN_FILLED, 'needs informative message') + _ASSERT(QLLS_FILLED, 'needs informative message') + _ASSERT(QICN_FILLED, 'needs informative message') + _ASSERT(QILS_FILLED, 'needs informative message') + _ASSERT(CLCN_FILLED, 'needs informative message') + _ASSERT(CLLS_FILLED, 'needs informative message') case (3) - ASSERT_(nn == 7) ! Q, QLCN, QLLS, QICN, QILS, QLIQ, QICE - ASSERT_(SPHU_FILLED) - ASSERT_(QLIQ_FILLED) - ASSERT_(QICE_FILLED) - ASSERT_(QLCN_FILLED) - ASSERT_(QLLS_FILLED) - ASSERT_(QICN_FILLED) - ASSERT_(QILS_FILLED) + _ASSERT(nn == 7, 'needs informative message') ! Q, QLCN, QLLS, QICN, QILS, QLIQ, QICE + _ASSERT(SPHU_FILLED, 'needs informative message') + _ASSERT(QLIQ_FILLED, 'needs informative message') + _ASSERT(QICE_FILLED, 'needs informative message') + _ASSERT(QLCN_FILLED, 'needs informative message') + _ASSERT(QLLS_FILLED, 'needs informative message') + _ASSERT(QICN_FILLED, 'needs informative message') + _ASSERT(QILS_FILLED, 'needs informative message') case (1) - ASSERT_(nn == 5) ! Q, QLCN, QLLS, QICN, QILS - ASSERT_(SPHU_FILLED) - ASSERT_(QLCN_FILLED) - ASSERT_(QLLS_FILLED) - ASSERT_(QICN_FILLED) - ASSERT_(QILS_FILLED) + _ASSERT(nn == 5, 'needs informative message') ! Q, QLCN, QLLS, QICN, QILS + _ASSERT(SPHU_FILLED, 'needs informative message') + _ASSERT(QLCN_FILLED, 'needs informative message') + _ASSERT(QLLS_FILLED, 'needs informative message') + _ASSERT(QICN_FILLED, 'needs informative message') + _ASSERT(QILS_FILLED, 'needs informative message') end select select case(FV_Atm(1)%flagstruct%nwat) @@ -1959,7 +1959,7 @@ subroutine FV_Run (STATE, CLOCK, GC, RC) endif enddo end select - ASSERT_(nn == FV_Atm(1)%ncnst) + _ASSERT(nn == FV_Atm(1)%ncnst, 'needs informative message') else do n=1,STATE%GRID%NQ nn=nn+1 @@ -1969,7 +1969,7 @@ subroutine FV_Run (STATE, CLOCK, GC, RC) state%vars%tracer(n)%content(:,:,:) = FV_Atm(1)%q(isc:iec,jsc:jec,1:npz,nn) endif enddo - ASSERT_(nn == FV_Atm(1)%ncnst) + _ASSERT(nn == FV_Atm(1)%ncnst, 'needs informative message') endif if (DEBUG) call debug_fv_state('After Dynamics Execution',STATE) @@ -3572,10 +3572,14 @@ subroutine fv_getUpdraftHelicity(uh25) FV_Atm(1)%bd%isd, FV_Atm(1)%bd%ied, FV_Atm(1)%bd%jsd, FV_Atm(1)%bd%jed, & FV_Atm(1)%npz, FV_Atm(1)%u, FV_Atm(1)%v, vort, & FV_Atm(1)%gridstruct%dx, FV_Atm(1)%gridstruct%dy, FV_Atm(1)%gridstruct%rarea) + +! call this with uh25_tmp which is of FVPRC call updraft_helicity(FV_Atm(1)%bd%isc, FV_Atm(1)%bd%iec, FV_Atm(1)%bd%jsc, FV_Atm(1)%bd%jec, FV_Atm(1)%ng, FV_Atm(1)%npz, & zvir, sphum, uh25_tmp, & FV_Atm(1)%w, vort, FV_Atm(1)%delz, FV_Atm(1)%q, & FV_Atm(1)%flagstruct%hydrostatic, FV_Atm(1)%pt, FV_Atm(1)%peln, FV_Atm(1)%phis, fms_grav, z_bot, z_top) + +! cast back to r4 uh25 = uh25_tmp end subroutine fv_getUpdraftHelicity diff --git a/GEOS_FV3_Utilities.F90 b/GEOS_FV3_Utilities.F90 new file mode 100644 index 0000000..8aecb06 --- /dev/null +++ b/GEOS_FV3_Utilities.F90 @@ -0,0 +1,244 @@ +module GEOS_FV3_UtilitiesMod + + implicit none + private + + public :: a2d2c + + + contains + + subroutine A2D2C(U,V,npz,getC) + +! Move A-Grid winds/tendencies oriented on lat/lon to the D-grid, or Optionally C-grid cubed-sphere orientation +! will return d-grid winds unless user asks for c-grid winds + + use mpp_domains_mod, only: mpp_update_domains, mpp_get_boundary, DGRID_NE + use mpp_parameter_mod, only: AGRID + use FV_StateMod, only: fv_atm + use fv_arrays_mod, only: REAL4, REAL8, FVPRC + use fv_mp_mod, only: is,js,ie,je,isd,jsd,ied,jed,ng + use sw_core_mod, only: d2a2c_vect + implicit none + + real, intent(INOUT) :: U(:,:,:) + real, intent(INOUT) :: V(:,:,:) + integer, intent( IN) :: npz + logical, intent( IN) :: getC + +!local variables + integer :: i,j,k, im2,jm2 + + real(REAL8) :: ud(is:ie,js:je+1,npz) + real(REAL8) :: vd(is:ie+1,js:je,npz) + + real(FVPRC) :: ut(isd:ied, jsd:jed) + real(FVPRC) :: vt(isd:ied, jsd:jed) + + real(REAL8) :: v3(is-1:ie+1,js-1:je+1,3) + real(REAL8) :: ue(is-1:ie+1,js :je+1,3) ! 3D winds at edges + real(REAL8) :: ve(is :ie+1,js-1:je+1,3) ! 3D winds at edges + real(REAL8), dimension(is:ie):: ut1, ut2, ut3 + real(REAL8), dimension(js:je):: vt1, vt2, vt3 + + real(FVPRC) :: uctemp(isd:ied+1,jsd:jed ,npz) + real(FVPRC) :: vctemp(isd:ied ,jsd:jed+1,npz) + real(FVPRC) :: utemp(isd:ied ,jsd:jed+1,npz) + real(FVPRC) :: vtemp(isd:ied+1,jsd:jed ,npz) + real(FVPRC) :: uatemp(isd:ied,jsd:jed,npz) + real(FVPRC) :: vatemp(isd:ied,jsd:jed,npz) +!#else +! real*8 :: uctemp(isd:ied+1,jsd:jed ,npz) +!! real*8 :: vctemp(isd:ied ,jsd:jed+1,npz) +!#endif + + real(FVPRC) :: wbuffer(js:je,npz) + real(FVPRC) :: sbuffer(is:ie,npz) + real(FVPRC) :: ebuffer(js:je,npz) + real(FVPRC) :: nbuffer(is:ie,npz) + integer :: npx, npy + integer :: STATUS + + npx = FV_Atm(1)%npx + npy = FV_Atm(1)%npy + + uatemp = 0.0d0 + vatemp = 0.0d0 + + uatemp(is:ie,js:je,:) = U + vatemp(is:ie,js:je,:) = V + + im2 = (npx-1)/2 + jm2 = (npy-1)/2 + +! Cubed-Sphere + call mpp_update_domains(uatemp, FV_Atm(1)%domain, complete=.false.) + call mpp_update_domains(vatemp, FV_Atm(1)%domain, complete=.true.) + do k=1, npz +! Compute 3D wind tendency on A grid + do j=js-1,je+1 + do i=is-1,ie+1 + v3(i,j,1) = uatemp(i,j,k)*fv_atm(1)%gridstruct%vlon(i,j,1) + vatemp(i,j,k)*fv_atm(1)%gridstruct%vlat(i,j,1) + v3(i,j,2) = uatemp(i,j,k)*fv_atm(1)%gridstruct%vlon(i,j,2) + vatemp(i,j,k)*fv_atm(1)%gridstruct%vlat(i,j,2) + v3(i,j,3) = uatemp(i,j,k)*fv_atm(1)%gridstruct%vlon(i,j,3) + vatemp(i,j,k)*fv_atm(1)%gridstruct%vlat(i,j,3) + enddo + enddo +! A --> D +! Interpolate to cell edges + do j=js,je+1 + do i=is-1,ie+1 + ue(i,j,1) = v3(i,j-1,1) + v3(i,j,1) + ue(i,j,2) = v3(i,j-1,2) + v3(i,j,2) + ue(i,j,3) = v3(i,j-1,3) + v3(i,j,3) + enddo + enddo + + do j=js-1,je+1 + do i=is,ie+1 + ve(i,j,1) = v3(i-1,j,1) + v3(i,j,1) + ve(i,j,2) = v3(i-1,j,2) + v3(i,j,2) + ve(i,j,3) = v3(i-1,j,3) + v3(i,j,3) + enddo + enddo +! --- E_W edges (for v-wind): + if ( is==1 ) then + i = 1 + do j=js,je + if ( j>jm2 ) then + vt1(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j-1,1)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,1) + vt2(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j-1,2)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,2) + vt3(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j-1,3)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,3) + else + vt1(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j+1,1)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,1) + vt2(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j+1,2)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,2) + vt3(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j+1,3)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,3) + endif + enddo + do j=js,je + ve(i,j,1) = vt1(j) + ve(i,j,2) = vt2(j) + ve(i,j,3) = vt3(j) + enddo + endif + if ( (ie+1)==npx ) then + i = npx + do j=js,je + if ( j>jm2 ) then + vt1(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j-1,1)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,1) + vt2(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j-1,2)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,2) + vt3(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j-1,3)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,3) + else + vt1(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j+1,1)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,1) + vt2(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j+1,2)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,2) + vt3(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j+1,3)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,3) + endif + enddo + do j=js,je + ve(i,j,1) = vt1(j) + ve(i,j,2) = vt2(j) + ve(i,j,3) = vt3(j) + enddo + endif +! N-S edges (for u-wind): + if ( js==1 ) then + j = 1 + do i=is,ie + if ( i>im2 ) then + ut1(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i-1,j,1)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,1) + ut2(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i-1,j,2)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,2) + ut3(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i-1,j,3)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,3) + else + ut1(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i+1,j,1)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,1) + ut2(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i+1,j,2)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,2) + ut3(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i+1,j,3)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,3) + endif + enddo + do i=is,ie + ue(i,j,1) = ut1(i) + ue(i,j,2) = ut2(i) + ue(i,j,3) = ut3(i) + enddo + endif + + if ( (je+1)==npy ) then + j = npy + do i=is,ie + if ( i>im2 ) then + ut1(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i-1,j,1)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,1) + ut2(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i-1,j,2)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,2) + ut3(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i-1,j,3)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,3) + else + ut1(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i+1,j,1)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,1) + ut2(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i+1,j,2)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,2) + ut3(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i+1,j,3)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,3) + endif + enddo + do i=is,ie + ue(i,j,1) = ut1(i) + ue(i,j,2) = ut2(i) + ue(i,j,3) = ut3(i) + enddo + endif +! Update: + do j=js,je+1 + do i=is,ie + ud(i,j,k) = 0.5*( ue(i,j,1)*fv_atm(1)%gridstruct%es(1,i,j,1) + & + ue(i,j,2)*fv_atm(1)%gridstruct%es(2,i,j,1) + & + ue(i,j,3)*fv_atm(1)%gridstruct%es(3,i,j,1) ) + enddo + enddo + do j=js,je + do i=is,ie+1 + vd(i,j,k) = 0.5*( ve(i,j,1)*fv_atm(1)%gridstruct%ew(1,i,j,2) + & + ve(i,j,2)*fv_atm(1)%gridstruct%ew(2,i,j,2) + & + ve(i,j,3)*fv_atm(1)%gridstruct%ew(3,i,j,2) ) + enddo + enddo + + enddo ! k-loop + + if (getC) then + + ! Now we have D-Grid winds, need to make call to d2a2c_vect + + utemp = 0.0d0 + vtemp = 0.0d0 + utemp(is:ie,js:je,:) = ud(is:ie,js:je,:) + vtemp(is:ie,js:je,:) = vd(is:ie,js:je,:) + + ! update shared edges + call mpp_get_boundary(utemp, vtemp, FV_Atm(1)%domain, & + wbuffery=wbuffer, ebuffery=ebuffer, & + sbufferx=sbuffer, nbufferx=nbuffer, & + gridtype=DGRID_NE, complete=.true. ) + do k=1,npz + do i=is,ie + utemp(i,je+1,k) = nbuffer(i,k) + enddo + do j=js,je + vtemp(ie+1,j,k) = ebuffer(j,k) + enddo + enddo + + call mpp_update_domains(utemp, vtemp, FV_Atm(1)%domain, gridtype=DGRID_NE, complete=.true.) + do k=1,npz + call d2a2c_vect(utemp(:,:,k), vtemp(:,:,k), & + uatemp(:,:,k), vatemp(:,:,k), & + uctemp(:,:,k), vctemp(:,:,k), ut, vt, .true., & + fv_atm(1)%gridstruct,fv_atm(1)%bd,npx,npy,.false.,0) + enddo + + U(:,:,:) = uctemp(is:ie,js:je,:) + V(:,:,:) = vctemp(is:ie,js:je,:) + + else + + ! return d-grid winds + u = ud(is:ie,js:je,:) + v = vd(is:ie,js:je,:) + + end if + + end subroutine A2D2C + +end module GEOS_FV3_UtilitiesMod diff --git a/GetWeights.F90 b/GetWeights.F90 index 8345ae2..c21f990 100644 --- a/GetWeights.F90 +++ b/GetWeights.F90 @@ -50,9 +50,7 @@ subroutine GetWeights(npx, npy, nlat, nlon, index, weight, id1, id2, jdc, l2c, & ee1, ee2, ff1, ff2, gg1, gg2, e1, e2, f1, f2, g1, g2, sublons, sublats, AmNodeRoot, WriteNetcdf) #include "MAPL_Generic.h" - use MAPL_BaseMod - use MAPL_GenericMod - use MAPL_ShmemMod + use MAPL use fv_grid_utils_mod, only : gnomonic_grids, cell_center2, mid_pt_sphere use fv_grid_tools_mod, only : mirror_grid use fv_grid_tools_mod, only : get_unit_vector @@ -368,7 +366,7 @@ subroutine GetWeights(npx, npy, nlat, nlon, index, weight, id1, id2, jdc, l2c, & if(.not. allocated(FV_Atm)) return call MAPL_SyncSharedMemory(rc=STATUS) - VERIFY_(STATUS) + _VERIFY(STATUS) allocate(e1(is:ie,js:je,3)) allocate(e2(is:ie,js:je,3)) allocate(f1(is:ie,js:je,3)) @@ -495,238 +493,3 @@ subroutine CreateCube2LatLonRotation(grid, center, ee1, ee2, ff1, ff2, gg1, gg2) end subroutine CreateCube2LatLonRotation end subroutine GetWeights - -! version of a2d3d routine from FV_StateMod - - subroutine A2D2C(U,V,npz,getC) - -! Move A-Grid winds/tendencies oriented on lat/lon to the D-grid, or Optionally C-grid cubed-sphere orientation -! will return d-grid winds unless user asks for c-grid winds - - use mpp_domains_mod, only: mpp_update_domains, mpp_get_boundary, DGRID_NE - use mpp_parameter_mod, only: AGRID - use FV_StateMod, only: fv_atm - use fv_arrays_mod, only: REAL4, REAL8, FVPRC - use fv_mp_mod, only: is,js,ie,je,isd,jsd,ied,jed,ng - use sw_core_mod, only: d2a2c_vect - implicit none - - real, intent(INOUT) :: U(:,:,:) - real, intent(INOUT) :: V(:,:,:) - integer, intent( IN) :: npz - logical, intent( IN) :: getC - -!local variables - integer :: i,j,k, im2,jm2 - - real(REAL8) :: ud(is:ie,js:je+1,npz) - real(REAL8) :: vd(is:ie+1,js:je,npz) - - real(FVPRC) :: ut(isd:ied, jsd:jed) - real(FVPRC) :: vt(isd:ied, jsd:jed) - - real(REAL8) :: v3(is-1:ie+1,js-1:je+1,3) - real(REAL8) :: ue(is-1:ie+1,js :je+1,3) ! 3D winds at edges - real(REAL8) :: ve(is :ie+1,js-1:je+1,3) ! 3D winds at edges - real(REAL8), dimension(is:ie):: ut1, ut2, ut3 - real(REAL8), dimension(js:je):: vt1, vt2, vt3 - - real(FVPRC) :: uctemp(isd:ied+1,jsd:jed ,npz) - real(FVPRC) :: vctemp(isd:ied ,jsd:jed+1,npz) - real(FVPRC) :: utemp(isd:ied ,jsd:jed+1,npz) - real(FVPRC) :: vtemp(isd:ied+1,jsd:jed ,npz) - real(FVPRC) :: uatemp(isd:ied,jsd:jed,npz) - real(FVPRC) :: vatemp(isd:ied,jsd:jed,npz) -!#else -! real*8 :: uctemp(isd:ied+1,jsd:jed ,npz) -!! real*8 :: vctemp(isd:ied ,jsd:jed+1,npz) -!#endif - - real(FVPRC) :: wbuffer(js:je,npz) - real(FVPRC) :: sbuffer(is:ie,npz) - real(FVPRC) :: ebuffer(js:je,npz) - real(FVPRC) :: nbuffer(is:ie,npz) - integer :: npx, npy - integer :: STATUS - - npx = FV_Atm(1)%npx - npy = FV_Atm(1)%npy - - uatemp = 0.0d0 - vatemp = 0.0d0 - - uatemp(is:ie,js:je,:) = U - vatemp(is:ie,js:je,:) = V - - im2 = (npx-1)/2 - jm2 = (npy-1)/2 - -! Cubed-Sphere - call mpp_update_domains(uatemp, FV_Atm(1)%domain, complete=.false.) - call mpp_update_domains(vatemp, FV_Atm(1)%domain, complete=.true.) - do k=1, npz -! Compute 3D wind tendency on A grid - do j=js-1,je+1 - do i=is-1,ie+1 - v3(i,j,1) = uatemp(i,j,k)*fv_atm(1)%gridstruct%vlon(i,j,1) + vatemp(i,j,k)*fv_atm(1)%gridstruct%vlat(i,j,1) - v3(i,j,2) = uatemp(i,j,k)*fv_atm(1)%gridstruct%vlon(i,j,2) + vatemp(i,j,k)*fv_atm(1)%gridstruct%vlat(i,j,2) - v3(i,j,3) = uatemp(i,j,k)*fv_atm(1)%gridstruct%vlon(i,j,3) + vatemp(i,j,k)*fv_atm(1)%gridstruct%vlat(i,j,3) - enddo - enddo -! A --> D -! Interpolate to cell edges - do j=js,je+1 - do i=is-1,ie+1 - ue(i,j,1) = v3(i,j-1,1) + v3(i,j,1) - ue(i,j,2) = v3(i,j-1,2) + v3(i,j,2) - ue(i,j,3) = v3(i,j-1,3) + v3(i,j,3) - enddo - enddo - - do j=js-1,je+1 - do i=is,ie+1 - ve(i,j,1) = v3(i-1,j,1) + v3(i,j,1) - ve(i,j,2) = v3(i-1,j,2) + v3(i,j,2) - ve(i,j,3) = v3(i-1,j,3) + v3(i,j,3) - enddo - enddo -! --- E_W edges (for v-wind): - if ( is==1 ) then - i = 1 - do j=js,je - if ( j>jm2 ) then - vt1(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j-1,1)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,1) - vt2(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j-1,2)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,2) - vt3(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j-1,3)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,3) - else - vt1(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j+1,1)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,1) - vt2(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j+1,2)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,2) - vt3(j) = fv_atm(1)%gridstruct%edge_vect_w(j)*ve(i,j+1,3)+(1.-fv_atm(1)%gridstruct%edge_vect_w(j))*ve(i,j,3) - endif - enddo - do j=js,je - ve(i,j,1) = vt1(j) - ve(i,j,2) = vt2(j) - ve(i,j,3) = vt3(j) - enddo - endif - if ( (ie+1)==npx ) then - i = npx - do j=js,je - if ( j>jm2 ) then - vt1(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j-1,1)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,1) - vt2(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j-1,2)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,2) - vt3(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j-1,3)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,3) - else - vt1(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j+1,1)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,1) - vt2(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j+1,2)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,2) - vt3(j) = fv_atm(1)%gridstruct%edge_vect_e(j)*ve(i,j+1,3)+(1.-fv_atm(1)%gridstruct%edge_vect_e(j))*ve(i,j,3) - endif - enddo - do j=js,je - ve(i,j,1) = vt1(j) - ve(i,j,2) = vt2(j) - ve(i,j,3) = vt3(j) - enddo - endif -! N-S edges (for u-wind): - if ( js==1 ) then - j = 1 - do i=is,ie - if ( i>im2 ) then - ut1(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i-1,j,1)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,1) - ut2(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i-1,j,2)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,2) - ut3(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i-1,j,3)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,3) - else - ut1(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i+1,j,1)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,1) - ut2(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i+1,j,2)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,2) - ut3(i) = fv_atm(1)%gridstruct%edge_vect_s(i)*ue(i+1,j,3)+(1.-fv_atm(1)%gridstruct%edge_vect_s(i))*ue(i,j,3) - endif - enddo - do i=is,ie - ue(i,j,1) = ut1(i) - ue(i,j,2) = ut2(i) - ue(i,j,3) = ut3(i) - enddo - endif - - if ( (je+1)==npy ) then - j = npy - do i=is,ie - if ( i>im2 ) then - ut1(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i-1,j,1)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,1) - ut2(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i-1,j,2)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,2) - ut3(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i-1,j,3)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,3) - else - ut1(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i+1,j,1)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,1) - ut2(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i+1,j,2)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,2) - ut3(i) = fv_atm(1)%gridstruct%edge_vect_n(i)*ue(i+1,j,3)+(1.-fv_atm(1)%gridstruct%edge_vect_n(i))*ue(i,j,3) - endif - enddo - do i=is,ie - ue(i,j,1) = ut1(i) - ue(i,j,2) = ut2(i) - ue(i,j,3) = ut3(i) - enddo - endif -! Update: - do j=js,je+1 - do i=is,ie - ud(i,j,k) = 0.5*( ue(i,j,1)*fv_atm(1)%gridstruct%es(1,i,j,1) + & - ue(i,j,2)*fv_atm(1)%gridstruct%es(2,i,j,1) + & - ue(i,j,3)*fv_atm(1)%gridstruct%es(3,i,j,1) ) - enddo - enddo - do j=js,je - do i=is,ie+1 - vd(i,j,k) = 0.5*( ve(i,j,1)*fv_atm(1)%gridstruct%ew(1,i,j,2) + & - ve(i,j,2)*fv_atm(1)%gridstruct%ew(2,i,j,2) + & - ve(i,j,3)*fv_atm(1)%gridstruct%ew(3,i,j,2) ) - enddo - enddo - - enddo ! k-loop - - if (getC) then - - ! Now we have D-Grid winds, need to make call to d2a2c_vect - - utemp = 0.0d0 - vtemp = 0.0d0 - utemp(is:ie,js:je,:) = ud(is:ie,js:je,:) - vtemp(is:ie,js:je,:) = vd(is:ie,js:je,:) - - ! update shared edges - call mpp_get_boundary(utemp, vtemp, FV_Atm(1)%domain, & - wbuffery=wbuffer, ebuffery=ebuffer, & - sbufferx=sbuffer, nbufferx=nbuffer, & - gridtype=DGRID_NE, complete=.true. ) - do k=1,npz - do i=is,ie - utemp(i,je+1,k) = nbuffer(i,k) - enddo - do j=js,je - vtemp(ie+1,j,k) = ebuffer(j,k) - enddo - enddo - - call mpp_update_domains(utemp, vtemp, FV_Atm(1)%domain, gridtype=DGRID_NE, complete=.true.) - do k=1,npz - call d2a2c_vect(utemp(:,:,k), vtemp(:,:,k), & - uatemp(:,:,k), vatemp(:,:,k), & - uctemp(:,:,k), vctemp(:,:,k), ut, vt, .true., & - fv_atm(1)%gridstruct,fv_atm(1)%bd,npx,npy,.false.,0) - enddo - - U(:,:,:) = uctemp(is:ie,js:je,:) - V(:,:,:) = vctemp(is:ie,js:je,:) - - else - - ! return d-grid winds - u = ud(is:ie,js:je,:) - v = vd(is:ie,js:je,:) - - end if - - end subroutine A2D2C diff --git a/LatLon2Cube.F90 b/LatLon2Cube.F90 index 9ff5438..035c8c2 100644 --- a/LatLon2Cube.F90 +++ b/LatLon2Cube.F90 @@ -1,7 +1,7 @@ subroutine latlon2cube(npx, npy, nlon, nlat, data_ll, data_cs) use ESMF - use MAPL_Mod, only : MAPL_UNDEF + use MAPL, only : MAPL_UNDEF use MAPL_IOMod, only : GETFILEUNIT, FREE_FILE use MAPL_ConstantsMod, only : pi=> MAPL_PI_R8 use fv_grid_utils_mod, only : gnomonic_grids, cell_center2 diff --git a/LatLonToCubeRegridder.F90 b/LatLonToCubeRegridder.F90 index 1f9508f..e002c45 100644 --- a/LatLonToCubeRegridder.F90 +++ b/LatLonToCubeRegridder.F90 @@ -5,10 +5,8 @@ #define _RETURN(A) if(present(rc)) rc=A; return module LatLonToCubeRegridderMod - use MAPL_AbstractRegridderMod + use MAPL use CubeLatLonTransformMod - use MAPL_GridSpecMod - use MAPL_RegridderSpecMod use, intrinsic :: iso_fortran_env, only: REAL32 use, intrinsic :: iso_fortran_env, only: REAL64 implicit none @@ -48,12 +46,7 @@ end function newLatLonToCubeRegridder subroutine initialize_subclass(this, unusable, rc) - use MAPL_KeywordEnforcerMod - use MAPL_BaseMod - use MAPL_CommsMod - use MAPL_AbstractGridFactoryMod - use MAPL_LatLonGridFactoryMod - use MAPL_GridManagerMod + use MAPL class (LatLonToCubeRegridder), intent(inout) :: this class (KeywordEnforcer), optional, intent(in) :: unusable integer, optional, intent(out) :: rc @@ -135,8 +128,7 @@ end subroutine regrid_scalar_2d_real32 subroutine regrid_scalar_3d_real32(this, q_in, q_out, rc) - use MAPL_CommsMod - use MAPL_BaseMod + use MAPL class (LatLonToCubeRegridder), intent(in) :: this real (kind=REAL32), intent(in) :: q_in(:,:,:) @@ -216,8 +208,7 @@ subroutine regrid_scalar_3d_real32(this, q_in, q_out, rc) end subroutine regrid_scalar_3d_real32 subroutine regrid_vector_3d_real32(this, u_in, v_in, u_out, v_out, rotate, rc) - use MAPL_CommsMod - use MAPL_BaseMod + use MAPL use, intrinsic :: iso_fortran_env, only: REAL32 class (LatLonToCubeRegridder), intent(in) :: this real(kind=REAL32), intent(in) :: u_in(:,:,:) @@ -303,8 +294,7 @@ subroutine transpose_regrid_scalar_2d_real32(this, q_in, q_out, rc) end subroutine transpose_regrid_scalar_2d_real32 subroutine transpose_regrid_scalar_3d_real32(this, q_in, q_out, rc) - use MAPL_CommsMod - use MAPL_BaseMod + use MAPL class (LatLonToCubeRegridder), intent(in) :: this real (kind=REAL32), intent(in) :: q_in(:,:,:) @@ -384,8 +374,7 @@ subroutine transpose_regrid_scalar_3d_real32(this, q_in, q_out, rc) end subroutine transpose_regrid_scalar_3d_real32 subroutine transpose_regrid_vector_3d_real32(this, u_in, v_in, u_out, v_out, rotate, rc) - use MAPL_CommsMod - use MAPL_BaseMod + use MAPL use, intrinsic :: iso_fortran_env, only: REAL32 class (LatLonToCubeRegridder), intent(in) :: this real(kind=REAL32), intent(in) :: u_in(:,:,:) diff --git a/StandAlone_AdvCore.F90 b/StandAlone_AdvCore.F90 index 77d6be1..4ee9473 100755 --- a/StandAlone_AdvCore.F90 +++ b/StandAlone_AdvCore.F90 @@ -3,7 +3,7 @@ #include "MAPL_Generic.h" program StandAlone_AdvCore - use MAPL_Mod + use MAPL use AdvCore_GridCompMod, only: SetServices use MPI diff --git a/StandAlone_DynAdvCore.F90 b/StandAlone_DynAdvCore.F90 index ee1b9f7..f2145c6 100755 --- a/StandAlone_DynAdvCore.F90 +++ b/StandAlone_DynAdvCore.F90 @@ -3,7 +3,7 @@ #include "MAPL_Generic.h" program StandAlone_DynAdvCore - use MAPL_Mod + use MAPL use StandAlone_DynAdvCore_GridCompMod, only: SetServices use MPI diff --git a/StandAlone_DynAdvCore_GridCompMod.F90 b/StandAlone_DynAdvCore_GridCompMod.F90 index d0b96b1..db57f4e 100755 --- a/StandAlone_DynAdvCore_GridCompMod.F90 +++ b/StandAlone_DynAdvCore_GridCompMod.F90 @@ -14,7 +14,7 @@ module StandAlone_DynAdvCore_GridCompMod ! !USES: use ESMF - use MAPL_Mod + use MAPL use AdvCore_GridCompMod, only : AdvCoreSetServices => SetServices use FVdycoreCubed_GridComp, only : DynCoreSetServices => SetServices diff --git a/StandAlone_FV3_Dycore.F90 b/StandAlone_FV3_Dycore.F90 index a98efc0..ffdd6de 100644 --- a/StandAlone_FV3_Dycore.F90 +++ b/StandAlone_FV3_Dycore.F90 @@ -4,10 +4,8 @@ #define _RC rc=status); _VERIFY(status program StandAlone_FV3_Dycore - use MAPL_Mod + use MAPL use FVdycoreCubed_GridComp, only: SetServices - use MAPL_CapOptionsMod - use MAPL_FlapCapOptionsMod implicit none !EOP diff --git a/c2c.F90 b/c2c.F90 index c42428c..ba90f86 100644 --- a/c2c.F90 +++ b/c2c.F90 @@ -1,7 +1,6 @@ -#define VERIFY_(A) if(MAPL_VRFY(A,Iam,__LINE__))call exit(-1) -#define ASSERT_(A) if(MAPL_ASRT(A,Iam,__LINE__))call exit(-1) - +#define I_AM_MAIN +#include "MAPL_Generic.h" #if defined(__INTEL_COMPILER) # define _FTELL ftelli8 @@ -14,8 +13,7 @@ program gmao_regrid use ESMF - use MAPL_Mod - use CubedSphereGridFactoryMod + use MAPL implicit none integer, parameter :: GridType_Unknown = 0 @@ -87,10 +85,12 @@ program gmao_regrid 2880, 1441], & ! G - 1/8 degree rectangular SHAPE = [2,8] ) - type(MAPL_NCIO) :: inNCIO,outNCIO + !type(MAPL_NCIO) :: inNCIO,outNCIO integer :: nDims, dimSizes(4),nSpatialDims type (ESMF_Grid) :: gridIn, gridOut + integer :: rc + ! Begin nargs = iargc() ! get command line argument info @@ -113,7 +113,7 @@ program gmao_regrid print *,'' print *,'ERROR: currently PARALLEL jobs not supported' print *,'' - ASSERT_(.false.) + _ASSERT(.false.,'needs informative message') end if if (MAPL_AM_I_Root(vm)) then @@ -129,10 +129,10 @@ program gmao_regrid ! determine grid type, and compute/guess im,im if (filetype ==0) then - InNCIO = MAPL_NCIOOpen(f_in,rc=status) - VERIFY_(STATUS) - call GetGridInfo(gi, filetype, ncinfo=InNCIO, rc=status) - VERIFY_(STATUS) + !InNCIO = MAPL_NCIOOpen(f_in,rc=status) + !VERIFY_(STATUS) + !call GetGridInfo(gi, filetype, ncinfo=InNCIO, rc=status) + !VERIFY_(STATUS) else call GetGridInfo(gi, filetype, rc=status) VERIFY_(STATUS) @@ -157,7 +157,7 @@ program gmao_regrid else if (str(1:1) == 'C') then if (str(2:2) /= ' ') then ic = ichar(str(2:2)) - ASSERT_(ic >= i0 .and. ic <= i9) ! Must be a number + _ASSERT(ic >= i0 .and. ic <= i9,'needs informative message') ! Must be a number read(str(2:),*) im gout%gridtype = GridType_CubedSphere gout%IM = im @@ -178,7 +178,7 @@ program gmao_regrid else if (str(1:1) == 'F') then if (str(2:2) /= ' ') then ic = ichar(str(2:2)) - ASSERT_(ic >= i0 .and. ic <= i9) ! Must be a number + _ASSERT(ic >= i0 .and. ic <= i9,'needs informative message') ! Must be a number read(str(2:),*) im gout%gridtype = GridType_CubedSphere gout%IM = im @@ -191,7 +191,7 @@ program gmao_regrid else if (str(1:1) == 'G') then if (str(2:2) /= ' ') then ic = ichar(str(2:2)) - ASSERT_(ic >= i0 .and. ic <= i9) ! Must be a number + _ASSERT(ic >= i0 .and. ic <= i9,'needs informative message') ! Must be a number read(str(2:),*) im gout%gridtype = GridType_CubedSphere gout%IM = im @@ -225,28 +225,28 @@ program gmao_regrid if (filetype ==0) then - call MAPL_NCIOChangeRes(InNCIO,OutNCIO,latSize=gout%jm,lonSize=gout%im,rc=status) - VERIFY_(STATUS) - call MAPL_NCIOSet(OutNCIO,filename=gout%filename) - call MAPL_NCIOCreateFile(OutNCIO) - do n=1,InNCIO%nVars - call MAPL_NCIOVarGetDims(InNCIO,InNCIO%vars(n)%name,nDims,dimSizes,nSpatialDims=nSpatialDims) + !call MAPL_NCIOChangeRes(InNCIO,OutNCIO,latSize=gout%jm,lonSize=gout%im,rc=status) + !VERIFY_(STATUS) + !call MAPL_NCIOSet(OutNCIO,filename=gout%filename) + !call MAPL_NCIOCreateFile(OutNCIO) + !do n=1,InNCIO%nVars + !call MAPL_NCIOVarGetDims(InNCIO,InNCIO%vars(n)%name,nDims,dimSizes,nSpatialDims=nSpatialDims) - if (nSpatialDims == 2) then - call MAPL_VarRead(InNCIO,InNCIO%vars(n)%name,var_in) - call regridder%regrid(var_in, var_out, rc=status) - VERIFY_(STATUS) - call MAPL_VarWrite(OutNCIO,InNCIO%vars(n)%name,var_out) - else if (nSpatialDims ==3) then - do i=1,dimSizes(3) - call MAPL_VarRead(InNCIO,InNCIO%vars(n)%name,var_in,lev=i) - call regridder%regrid(var_in, var_out, rc=status) - VERIFY_(STATUS) - call MAPL_VarWrite(OutNCIO,InNCIO%vars(n)%name,var_out,lev=i) - end do - end if - enddo - call MAPL_NCIOClose(OutNCIO) + !if (nSpatialDims == 2) then + !call MAPL_VarRead(InNCIO,InNCIO%vars(n)%name,var_in) + !call regridder%regrid(var_in, var_out, rc=status) + !VERIFY_(STATUS) + !call MAPL_VarWrite(OutNCIO,InNCIO%vars(n)%name,var_out) + !else if (nSpatialDims ==3) then + !do i=1,dimSizes(3) + !call MAPL_VarRead(InNCIO,InNCIO%vars(n)%name,var_in,lev=i) + !call regridder%regrid(var_in, var_out, rc=status) + !VERIFY_(STATUS) + !call MAPL_VarWrite(OutNCIO,InNCIO%vars(n)%name,var_out,lev=i) + !end do + !end if + !enddo + !call MAPL_NCIOClose(OutNCIO) else ! open files UNIT_R = GetFile(gi%filename, rc=status) @@ -275,7 +275,7 @@ program gmao_regrid RecEnd = _FTELL(unit_r) backspace(unit_r) RecStart = _FTELL(unit_r) - ASSERT_(4*((LM+1)+2) == RecEnd-RecStart) + _ASSERT(4*((LM+1)+2) == RecEnd-RecStart,'needs informative message') print *,'WARNING: encoutered shorter record, assuming PREF' read (unit_r) pref write(unit_w) pref @@ -296,18 +296,15 @@ program gmao_regrid call ESMF_Finalize (RC=status) VERIFY_(STATUS) -#undef VERIFY_ -#undef ASSERT_ - -#include "MAPL_Generic.h" - contains ! ================================================================ +#undef I_AM_MAIN + subroutine GuessFileType(filename, filetype, rc) use ESMF - use MAPL_Mod + use MAPL implicit none @@ -378,16 +375,16 @@ subroutine GuessFileType(filename, filetype, rc) end subroutine GuessFileType - subroutine GetGridInfo(gi, filetype, ncinfo, rc) + !subroutine GetGridInfo(gi, filetype, ncinfo, rc) + subroutine GetGridInfo(gi, filetype, rc) use ESMF - use MAPL_Mod - use MAPL_IOMod + use MAPL implicit none type(Regrid_GridInfo) :: gi integer :: filetype - type(MAPL_NCIO), optional, intent(in) :: ncinfo + !type(MAPL_NCIO), optional, intent(in) :: ncinfo integer, optional, intent(OUT):: RC integer :: i6, im, jm @@ -396,23 +393,23 @@ subroutine GetGridInfo(gi, filetype, ncinfo, rc) gi%gridtype = GridType_Unknown if (filetype == 0) then - gi%im=-1 - gi%jm=-1 - do i=1,ncinfo%ndims - if ( trim(ncinfo%dims(i)%name) == "lon" ) then - gi%im = ncinfo%dims(i)%len - else if (trim(ncinfo%dims(i)%name) == "lat" ) then - gi%jm = ncinfo%dims(i)%len - end if - enddo - ASSERT_(gi%im /= -1) - ASSERT_(gi%jm /= -1) - if (gi%jm == gi%im*6) then - gi%gridtype = GridType_CubedSphere - else - gi%gridtype = GridType_LatLon - end if - RETURN_(ESMF_SUCCESS) + !gi%im=-1 + !gi%jm=-1 + !do i=1,ncinfo%ndims + !if ( trim(ncinfo%dims(i)%name) == "lon" ) then + !gi%im = ncinfo%dims(i)%len + !else if (trim(ncinfo%dims(i)%name) == "lat" ) then + !gi%jm = ncinfo%dims(i)%len + !end if + !enddo + !_ASSERT(gi%im /= -1,'needs informative message') + !_ASSERT(gi%jm /= -1,'needs informative message') + !if (gi%jm == gi%im*6) then + !gi%gridtype = GridType_CubedSphere + !else + !gi%gridtype = GridType_LatLon + !end if + !RETURN_(ESMF_SUCCESS) end if if (filetype == 6) then diff --git a/c2c_rst.F90 b/c2c_rst.F90 index faca643..38bea26 100644 --- a/c2c_rst.F90 +++ b/c2c_rst.F90 @@ -1,6 +1,6 @@ #define VERIFY_(A) if(MAPL_VRFY(A,Iam,__LINE__))call exit(-1) -#define ASSERT_(A) if(MAPL_ASRT(A,Iam,__LINE__))call exit(-1) +#define _ASSERT(A) if(MAPL_ASRT(A,Iam,__LINE__))call exit(-1,'needs informative message') #if defined(__INTEL_COMPILER) @@ -14,7 +14,7 @@ program gmao_regrid use ESMF - use MAPL_Mod + use MAPL use, intrinsic :: iso_fortran_env implicit none @@ -118,7 +118,7 @@ program gmao_regrid print *,'' print *,'ERROR: currently PARALLEL jobs not supported' print *,'' - ASSERT_(.false.) + _ASSERT(.false.,'needs informative message') end if if (MAPL_AM_I_Root(vm)) then @@ -162,7 +162,7 @@ program gmao_regrid else if (str(1:1) == 'C') then if (str(2:2) /= ' ') then ic = ichar(str(2:2)) - ASSERT_(ic >= i0 .and. ic <= i9) ! Must be a number + _ASSERT(ic >= i0 .and. ic <= i9,'needs informative message') ! Must be a number read(str(2:),*) im gout%gridtype = GridType_CubedSphere gout%IM = im @@ -186,7 +186,7 @@ program gmao_regrid changeResolution = gi%im /= gout%im .or. gi%jm /= gout%jm if (gout%gridtype /= GridType_CubedSphere .and. gi%gridtype /= GridType_CubedSphere) then - ASSERT_(.false.) + _ASSERT(.false.,'needs informative message') end if call MAPL_GenGridName(gi%im, gi%jm, gridname=gridnamef, geos_style=.false.) call MAPL_GenGridName(gout%im, gout%jm, gridname=gridnameo, geos_style=.false.) @@ -348,7 +348,7 @@ program gmao_regrid RecEnd = _FTELL(unit_r) backspace(unit_r) RecStart = _FTELL(unit_r) - ASSERT_(4*((LM+1)+2) == RecEnd-RecStart) + _ASSERT(4*((LM+1)+2) == RecEnd-RecStart,'needs informative message') print *,'WARNING: encoutered shorter record, assuming PREF' read (unit_r) pref write(unit_w) pref @@ -374,7 +374,7 @@ program gmao_regrid VERIFY_(STATUS) #undef VERIFY_ -#undef ASSERT_ +#undef _ASSERT #include "MAPL_Generic.h" @@ -384,7 +384,7 @@ program gmao_regrid subroutine GuessFileType(filename, filetype, rc) use ESMF - use MAPL_Mod + use MAPL implicit none @@ -457,8 +457,7 @@ end subroutine GuessFileType subroutine GetGridInfo(gi, filetype, filename, ncinfo, rc) use ESMF - use MAPL_Mod - use MAPL_IOMod + use MAPL implicit none @@ -486,8 +485,8 @@ subroutine GetGridInfo(gi, filetype, filename, ncinfo, rc) gi%jm = ncinfo%dims(i)%len end if enddo - ASSERT_(gi%im /= -1) - ASSERT_(gi%jm /= -1) + _ASSERT(gi%im /= -1,'needs informative message') + _ASSERT(gi%jm /= -1,'needs informative message') if (gi%jm == gi%im*6) then gi%gridtype = GridType_CubedSphere else diff --git a/c2l_CFIO_offline.F90 b/c2l_CFIO_offline.F90 index 94a94fe..625978f 100644 --- a/c2l_CFIO_offline.F90 +++ b/c2l_CFIO_offline.F90 @@ -20,7 +20,7 @@ Program c2l_CFIO_offline use ESMF - use MAPL_Mod + use MAPL implicit NONE diff --git a/interp_restarts.F90 b/interp_restarts.F90 index 62d2243..a8b74de 100755 --- a/interp_restarts.F90 +++ b/interp_restarts.F90 @@ -1,7 +1,6 @@ #define VERIFY_(A) if((A)/=0) then; PRINT *, 'Interp_restarts.x', __LINE__; call MPI_Abort(A); endif program interp_restarts -#define DEALLOCGLOB_(A) if(associated(A)) then;A=0;call MAPL_DeAllocNodeArray(A,rc=STATUS);if(STATUS==MAPL_NoShm) deallocate(A,stat=STATUS);NULLIFY(A);endif !--------------------------------------------------------------------! ! purpose: driver for interpolation of GEOS FV and Moist restarts ! @@ -22,12 +21,10 @@ program interp_restarts ! use fv_eta_mod, only: set_eta use m_set_eta, only: set_eta use memutils_mod, only: print_memuse_stats - use MAPL_IOMod - use MAPL_ShmemMod - use MAPL_ConstantsMod + use MAPL + use gFTL_StringVector + use gFTL_StringIntegerMap use rs_scaleMod - use pFIO_StringVectorMod - use pFIO_StringIntegerMapMod implicit none @@ -69,22 +66,36 @@ program interp_restarts integer :: nmoist logical :: isBinFV, isBinMoist integer :: ftype - type(MAPL_NCIO) :: ncio,ncioOut + type(Netcdf4_Fileformatter) :: InFmt,OutFmt + type(FileMetadata), allocatable :: InCfg(:),OutCfg(:) integer :: nVars,imc,jmc,lonid,latid,levid,edgeid character(62) :: vname - type(StringVector) :: moist_variables + type(StringVector) :: moist_variables,all_moist_vars + type(StringVectorIterator) :: siter type(StringIntegerMap) :: tracer_names + ! bma added character(len=128) :: moist_order(9) = (/"Q ","QLLS","QLCN","CLLS","CLCN","QILS","QICN","NCPL","NCPI"/) - integer :: p_split, npx, npy, npz,ivar,lcnt_var,iq0 + integer :: p_split, npx, npy, npz, ivar, lcnt_var, iq0 integer :: n_args,n_files,nlevs,nedges,ifile,nlev,n_output,nfv_vars character(len=ESMF_MAXPATHLEN), allocatable :: extra_files(:),extra_output(:) type(fv_rst), pointer :: rst_files(:) => null() type(ArrDescr) :: ArrDes integer :: info logical :: amWriter - integer :: isl,iel,jsl,jel,npes_x,npes_y,n_writers,n_readers,ndims,dimsize(3) + integer :: isl,iel,jsl,jel,npes_x,npes_y,n_writers,n_readers type(ESMF_Grid) :: grid logical :: scale_rst + type(StringVariableMap), pointer :: vars + type(StringVariableMapIterator) :: iter + character(len=:), pointer :: var_name + type(StringVariableMap), pointer :: variables + type(Variable), pointer :: myVariable + type(StringVariableMapIterator) :: var_iter + type(StringVector), pointer :: var_dimensions + character(len=:), pointer :: dname + integer :: dim1,ndims + type(CubedSphereGridFactory) :: csfactory + real, allocatable :: schmidt_parameters(:) real(FVPRC), allocatable :: q_wat(:,:,:,:) @@ -158,6 +169,14 @@ program interp_restarts else write(*,*)'bad argument to scalers, will scale by default' end if + case('-stretched_grid') + allocate(schmidt_parameters(3)) + call get_command_argument(i+1,astr) + read(astr,*)schmidt_parameters(1) + call get_command_argument(i+2,astr) + read(astr,*)schmidt_parameters(2) + call get_command_argument(i+3,astr) + read(astr,*)schmidt_parameters(3) end select end do @@ -174,6 +193,12 @@ program interp_restarts if (ihydro == 0) FV_Atm(1)%flagstruct%hydrostatic = .false. FV_Atm(1)%flagstruct%Make_NH = .false. if (.not. FV_Atm(1)%flagstruct%hydrostatic) FV_Atm(1)%flagstruct%Make_NH = .true. + if (allocated(schmidt_parameters)) then + FV_Atm(1)%flagstruct%do_schmidt = .true. + FV_Atm(1)%flagstruct%target_lon=schmidt_parameters(1) + FV_Atm(1)%flagstruct%target_lat=schmidt_parameters(2) + FV_Atm(1)%flagstruct%stretch_fac=schmidt_parameters(3) + end if if (n_files > 0) allocate(rst_files(n_files)) @@ -199,6 +224,16 @@ program interp_restarts isBinFV = .true. isBinMoist = .true. +! Determine Total Number of Tracers (MOIST, GOCART, PCHEM, ANA) +! ------------------------------------------------------------- + nmoist = 0 + isBinFV = .true. + isBinMoist = .true. + + call print_memuse_stats('interp_restarts: rs_count') + call mpp_broadcast(nmoist, mpp_root_pe()) + call mpp_broadcast(isBinMoist, mpp_root_pe()) + if (is_master()) print*, 'HYDROSTATIC : ', FV_Atm(1)%flagstruct%hydrostatic if (is_master()) print*, 'Make_NH : ', FV_Atm(1)%flagstruct%Make_NH if (is_master()) print*, 'Tracers : ', FV_Atm(1)%ncnst @@ -221,10 +256,15 @@ program interp_restarts km=header(3) close(IUNIT) else - ncio = MAPL_NCIOOpen("fvcore_internal_restart_in") - call MAPL_NCIOGetDimSizes(NCIO,lon=im,lat=jm,lev=km) - call MAPL_NCIOClose(NCIO,destroy=.true.) - end if + call InFmt%open("fvcore_internal_restart_in",pFIO_READ,rc=status) + allocate(InCfg(1)) + InCfg(1) = InFmt%read() + im = InCfg(1)%get_dimension('lon') + jm = InCfg(1)%get_dimension('lat') + km = InCfg(1)%get_dimension('lev') + call InFmt%close() + deallocate(InCfg) + end if else call mpp_error(FATAL, 'ABORT: fvcore_internal_restart_in does not exist') endif @@ -233,30 +273,39 @@ program interp_restarts call MAPL_NCIOGetFileType("moist_internal_restart_in",ftype) if (ftype == 0) then isBinMoist = .false. - NCIO = MAPL_NCIOOpen("moist_internal_restart_in") - call MAPL_NCIOGetDimSizes(ncio,slices=nmoist,nVars=nVars) + call InFmt%open("moist_internal_restart_in",pFIO_READ,rc=status) + allocate(InCfg(1)) + InCfg(1) = InFmt%read() + call MAPL_IOCountLevels(InCfg(1),nmoist) + all_moist_vars = MAPL_IOGetNonDimVars(InCfg(1),rc=status) + siter = all_moist_vars%begin() + variables => InCfg(1)%get_variables() lcnt_var=2 - do ivar=1,nVars - call MAPL_NCIOGetVarName(NCIO,ivar,vname) - call MAPL_NCIOVarGetDims(NCIO,vname,ndims,dimSize) + do while (siter /= all_moist_vars%end()) + var_name => siter%get() + myVariable => variables%at(var_name) + var_dimensions => myVariable%get_dimensions() + ndims = var_dimensions%size() if (ndims==2) nmoist=nmoist-1 if (ndims==3) then - call moist_variables%push_back(trim(vname)) - if (trim(vname)=='Q') then + call moist_variables%push_back(trim(var_name)) + if (trim(var_name)=='Q') then iq0=1 - call tracer_names%insert(trim(vname),iq0) + call tracer_names%insert(trim(var_name),iq0) else iq0=lcnt_var - call tracer_names%insert(trim(vname),iq0) + call tracer_names%insert(trim(var_name),iq0) lcnt_var=lcnt_var+1 end if end if - enddo - call MAPL_NCIOClose(ncio,destroy=.true.) + call siter%next() + end do + call InFmt%close() + deallocate(InCfg) else call rs_count( "moist_internal_restart_in",nmoist ) - if(mod(nmoist,km)/=0) then - call mpp_error(FATAL, 'ABORT: '//'binary moist restart must have only 3D variabels') + if (mod(nmoist,km)/=0) then + call mpp_error(FATAL, 'ABORT: '//'binary moist restart must have only 3D variabels') end if nVars = nmoist/km do ivar=1,nVars @@ -270,6 +319,8 @@ program interp_restarts else call mpp_error(FATAL, 'ABORT: moist_internal_restart_in does not exist') endif + + call print_memuse_stats('interp_restarts: rs_count') call mpp_broadcast(nmoist, mpp_root_pe()) call mpp_broadcast(isBinMoist, mpp_root_pe()) @@ -307,27 +358,44 @@ program interp_restarts if (ftype ==0) then rst_files(i)%isBin=.false. rst_files(i)%file_name=trim(extra_files(i)) - NCIO = MAPL_NCIOOpen(trim(extra_files(i))) - allocate(rst_files(i)%vars(size(ncio%vars))) - edgeid=-1 - levid=-1 - call MAPL_NCIOGetDimSizes(ncio,lev=nLevs,edges=nEdges,levid=levid,edgeid=edgeid) - do n=1,size(ncio%vars) - rst_files(i)%vars(n)%name=trim(ncio%vars(n)%name) - if (ncio%vars(n)%ndims==3) then - do j=1,size(ncio%vars(n)%dimids) - if (ncio%vars(n)%dimids(j) == edgeid) then - rst_files(i)%vars(n)%nlev= npz+1 - else if (ncio%vars(n)%dimids(j) == levid) then - rst_files(i)%vars(n)%nlev= npz + + call InFmt%open(trim(extra_files(i)),pFIO_READ,rc=status) + allocate(InCfg(1)) + InCfg(1) = InFmt%read() + call MAPL_IOCountNonDimVars(InCfg(1),nVars) + + allocate(rst_files(i)%vars(nVars)) + + variables => InCfg(1)%get_variables() + var_iter = variables%begin() + n=0 + do while (var_iter /= variables%end()) + + var_name => var_iter%key() + myVariable => var_iter%value() + if (.not.InCfg(1)%is_coordinate_variable(var_name)) then + n=n+1 + var_dimensions => myVariable%get_dimensions() + ndims = var_dimensions%size() + rst_files(i)%vars(n)%name=trim(var_name) + if (ndims ==2) then + rst_files(i)%vars(n)%nlev=1 + else if (ndims==3) then + dname => myVariable%get_ith_dimension(3) + dim1=InCfg(1)%get_dimension(dname) + if (dim1 == km) then + rst_files(i)%vars(n)%nlev=npz + else if (dim1 == km+1) then + rst_files(i)%vars(n)%nlev=npz+1 end if - enddo - else if (ncio%vars(n)%ndims==2) then - rst_files(i)%vars(n)%nlev=1 + end if end if + call var_iter%next() enddo + rst_files(i)%have_descriptor=.true. - call MAPL_NCIOClose(NCIO,destroy=.true.) + call InFmt%close() + deallocate(InCfg) else call rs_count(trim(extra_files(i)),nlevs) if (mod(nlevs,km) /= 0) then @@ -354,19 +422,35 @@ program interp_restarts call print_memuse_stats('interp_restarts: begining get_external_ic') + npes_x=fv_atm(1)%layout(1) + npes_y=fv_atm(1)%layout(2) + is = FV_Atm(1)%bd%isc + ie = FV_Atm(1)%bd%iec + js = FV_Atm(1)%bd%jsc + je = FV_Atm(1)%bd%jec + isl=is + iel=ie + jsl=(npx-1)*(tile-1)+js + jel=(npx-1)*(tile-1)+je + + call ArrDescrInit(Arrdes,MPI_COMM_WORLD,npx-1,(npx-1)*6,npz,npes_x,npes_y*6,n_readers,n_writers,isl,iel,jsl,jel,rc=status) + call ArrDescrSet(arrdes,offset=0_MPI_OFFSET_KIND) + if (allocated(schmidt_parameters)) then + csfactory = CubedSphereGridFactory(im_world=npx-1,lm=npz,nx=npes_x,ny=npes_y,stretch_factor=schmidt_parameters(3), & + target_lon=schmidt_parameters(1),target_lat=schmidt_parameters(2)) + else + csfactory = CubedSphereGridFactory(im_world=npx-1,lm=npz,nx=npes_x,ny=npes_y) + end if + grid = grid_manager%make_grid(csfactory,rc=status) + FV_Atm(1)%flagstruct%Make_NH = .false. ! Do this after rescaling if (jm == 6*im) then - call get_geos_ic( FV_Atm, rst_files, .true.) + call get_geos_ic( FV_Atm, rst_files, .true., grid) else - call get_geos_ic( FV_Atm, rst_files, .false.) + call get_geos_ic( FV_Atm, rst_files, .false., grid) endif FV_Atm(1)%flagstruct%Make_NH = .true. ! Reset this for later - is = FV_Atm(1)%bd%isc - ie = FV_Atm(1)%bd%iec - js = FV_Atm(1)%bd%jsc - je = FV_Atm(1)%bd%jec - if (scale_rst) then call scale_drymass(fv_atm,tracer_names,rc=status) VERIFY_(status) @@ -393,19 +477,7 @@ program interp_restarts pt_local(is:ie,js:je,k) = FV_Atm(1)%pt(is:ie,js:je,k)/FV_Atm(1)%pkz(is:ie,js:je,k) enddo - - - npes_x=fv_atm(1)%layout(1) - npes_y=fv_atm(1)%layout(2) - isl=is - iel=ie - jsl=(npx-1)*(tile-1)+js - jel=(npx-1)*(tile-1)+je - - call ArrDescrInit(Arrdes,MPI_COMM_WORLD,npx-1,(npx-1)*6,npes_x,npes_y*6,n_readers,n_writers,isl,iel,jsl,jel,rc=status) - call ArrDescrSet(arrdes,offset=0_MPI_OFFSET_KIND) amWriter = arrdes%writers_comm/=MPI_COMM_NULL - grid = simple_cs_grid_creator(npx-1,(npx-1)*6,npes_x,npes_y*6,status) call MPI_Info_create(info,status) call print_memuse_stats('interp_restarts: going to write restarts') @@ -413,189 +485,206 @@ program interp_restarts ! write fvcore_internal_rst if( file_exist("fvcore_internal_restart_in") ) then - write(fname1, "('fvcore_internal_rst_c',i4.4,'_',i3.3,'L')") npx-1,npz - if (is_master()) print*, 'Writing : ', TRIM(fname1) - if (isBinFV) then + write(fname1, "('fvcore_internal_rst_c',i4.4,'_',i3.3,'L')") npx-1,npz + if (is_master()) print*, 'Writing : ', TRIM(fname1) + if (isBinFV) then - open(IUNIT,file='fvcore_internal_restart_in' ,access='sequential',form='unformatted',status='old') - !if (AmWriter) open(OUNIT,file=TRIM(fname1),access='sequential',form='unformatted') - if (n_writers==1) then - OUNIT=getfile(TRIM(fname1),form='unformatted',rc=status) - VERIFY_(status) - else - if (AmWriter) then - call MPI_FILE_OPEN(arrdes%writers_comm, fname1, MPI_MODE_WRONLY+MPI_MODE_CREATE, & - info, OUNIT, STATUS) - VERIFY_(STATUS) - end if + open(IUNIT,file='fvcore_internal_restart_in' ,access='sequential',form='unformatted',status='old') + !if (AmWriter) open(OUNIT,file=TRIM(fname1),access='sequential',form='unformatted') + if (n_writers==1) then + OUNIT=getfile(TRIM(fname1),form='unformatted',rc=status) + VERIFY_(status) + else + if (AmWriter) then + call MPI_FILE_OPEN(arrdes%writers_comm, fname1, MPI_MODE_WRONLY+MPI_MODE_CREATE, & + info, OUNIT, STATUS) + VERIFY_(STATUS) end if + end if - ! Headers - read (IUNIT, IOSTAT=status) header - if(n_writers > 1) then - call Write_Parallel(HEADER, OUNIT, ARRDES=ARRDES, RC=status) - VERIFY_(STATUS) - else - if (amwriter) write(OUNIT) header - endif - if (is_master()) print*, header + ! Headers + read (IUNIT, IOSTAT=status) header + if(n_writers > 1) then + call Write_Parallel(HEADER, OUNIT, ARRDES=ARRDES, RC=status) + VERIFY_(STATUS) + else + if (amwriter) write(OUNIT) header + endif + if (is_master()) print*, header - read (IUNIT, IOSTAT=status) header(1:5) - if (is_master()) print*, header(1:5) - header(1) = (npx-1) - header(2) = (npy-1)*6 - header(3) = npz + read (IUNIT, IOSTAT=status) header(1:5) + if (is_master()) print*, header(1:5) + header(1) = (npx-1) + header(2) = (npy-1)*6 + header(3) = npz + + if(n_writers > 1) then + call Write_Parallel(HEADER(1:5), OUNIT, ARRDES=ARRDES, RC=status) + VERIFY_(STATUS) + else + if (amwriter) write(OUNIT) header(1:5) + endif - if(n_writers > 1) then - call Write_Parallel(HEADER(1:5), OUNIT, ARRDES=ARRDES, RC=status) - VERIFY_(STATUS) - else - if (amwriter) write(OUNIT) header(1:5) - endif + if (is_master()) print*, header(1:5) + close(IUNIT) - if (is_master()) print*, header(1:5) - close(IUNIT) + else - else + call InFmt%open("fvcore_internal_restart_in",pFIO_READ,rc=status) + allocate(InCfg(1),OutCfg(1)) + InCfg(1)=InFmt%read(rc=status) + call MAPL_IOCountNonDimVars(InCfg(1),nVars) - ncio = MAPL_NCIOOpen("fvcore_internal_restart_in",rc=status) - call MAPL_NCIOGetDimSizes(ncio,nVars=nVars) - if (AmWriter) then - imc = npx-1 - jmc = imc*6 - call MAPL_NCIOChangeRes(ncio,ncioOut,lonSize=imc,latSize=jmc,levSize=npz) - if (.not.fv_atm(1)%flagstruct%hydrostatic) then - nfv_vars = 9 - else - nfv_vars = 7 - end if - call MAPL_NCIOSet(ncioOut,filename=fname1,nVars=nfv_vars,overwriteVars=.true.) - !start creating file. Can not simply use input because if it is lat-lon will not have dz or w - call MAPL_NCIOGetDimSizes(ncioOut,lonid=lonid,latid=latid,levid=levid,edgeid=edgeid) - call MAPL_NCIOAddVar(ncioOut,"AK",(/edgeid/),6,units="Pa",long_name="hybrid_sigma_pressure_a",rc=status) - VERIFY_(status) - call MAPL_NCIOAddVar(ncioOut,"BK",(/edgeid/),6,units="Pa",long_name="hybrid_sigma_pressure_b",rc=status) - VERIFY_(status) - call MAPL_NCIOAddVar(ncioOut,"U",(/lonid,latid,levid/),6,units="m s-1",long_name="eastward_wind",rc=status) - VERIFY_(status) - call MAPL_NCIOAddVar(ncioOut,"V",(/lonid,latid,levid/),6,units="m s-1",long_name="northward_wind",rc=status) - VERIFY_(status) - call MAPL_NCIOAddVar(ncioOut,"PT",(/lonid,latid,levid/),6,units="K Pa$^{-\kappa}$",long_name="scaled_potential_temperature",rc=status) - VERIFY_(status) - call MAPL_NCIOAddVar(ncioOut,"PE",(/lonid,latid,edgeid/),6,units="Pa",long_name="air_pressure",rc=status) - VERIFY_(status) - call MAPL_NCIOAddVar(ncioOut,"PKZ",(/lonid,latid,levid/),6,units="Pa$^\kappa$",long_name="pressure_to_kappa",rc=status) - VERIFY_(status) - - ! if dz and w were not in the original file add them - ! they need to be in there for the restart - - if (.not.fv_atm(1)%flagstruct%hydrostatic) then - call MAPL_NCIOAddVar(ncioOut,"DZ",(/lonid,latid,levid/),6,units="m",long_name="height_thickness",rc=status) - VERIFY_(status) - call MAPL_NCIOAddVar(ncioOut,"W",(/lonid,latid,levid/),6,units="m s-1",long_name="vertical_velocity",rc=status) - VERIFY_(status) - endif - call MAPL_NCIOCreateFile(ncioOut,comm=arrdes%writers_comm,info=info,rc=status) - VERIFY_(status) - end if + if (AmWriter) then + imc = npx-1 + jmc = imc*6 + call MAPL_IOChangeRes(InCfg(1),OutCfg(1),(/'lon ','lat ','lev ','edge'/),(/imc,jmc,npz,npz+1/),rc=status) + + ! if dz and w were not in the original file add them + ! they need to be in there for the restart + if (.not.fv_atm(1)%flagstruct%hydrostatic) then + ! fix thic + !call MAPL_NCIOAddVar(ncioOut,"DZ",(/lonid,latid,levid/),6,units="m",long_name="height_thickness",rc=status) + !call MAPL_NCIOAddVar(ncioOut,"W",(/lonid,latid,levid/),6,units="m s-1",long_name="vertical_velocity",rc=status) + endif + call OutFmt%create_par(fname1,comm=arrdes%writers_comm,info=info,rc=status) + call OutFmt%write(OutCfg(1),rc=status) end if + end if + ! AK and BK - allocate ( r8_akbk(npz+1) ) - r8_akbk = FV_Atm(1)%ak - if (isBinFV) then - if (n_writers == 1) then - if (AmWriter) write(OUNIT) r8_akbk - else - call write_parallel(r8_akbk,ounit,arrdes=arrdes,rc=status) - VERIFY_(status) - end if - else - if (AmWriter) call MAPL_VarWrite(ncioOut,"AK",r8_akbk) + allocate ( r8_akbk(npz+1) ) + r8_akbk = FV_Atm(1)%ak + if (isBinFV) then + if (n_writers == 1) then + if (AmWriter) write(OUNIT) r8_akbk + else + call write_parallel(r8_akbk,ounit,arrdes=arrdes,rc=status) + VERIFY_(status) end if - r8_akbk = FV_Atm(1)%bk - if (isBinFV) then - if (n_writers == 1) then - if (AmWriter) write(OUNIT) r8_akbk - else - call write_parallel(r8_akbk,ounit,arrdes=arrdes,rc=status) - VERIFY_(status) - end if - else - if (AmWriter) call MAPL_VarWrite(ncioOut,"BK",r8_akbk) + else + write(*,*)'bma writing ak' + if (AmWriter) call MAPL_VarWrite(OutFmt,"AK",r8_akbk) + end if + r8_akbk = FV_Atm(1)%bk + if (isBinFV) then + if (n_writers == 1) then + if (AmWriter) write(OUNIT) r8_akbk + else + call write_parallel(r8_akbk,ounit,arrdes=arrdes,rc=status) + VERIFY_(status) end if - deallocate ( r8_akbk ) + else + write(*,*)'bma writing bk' + if (AmWriter) call MAPL_VarWrite(OutFmt,"BK",r8_akbk) + end if + deallocate ( r8_akbk ) - allocate(r8_local(is:ie,js:je,npz+1)) + allocate(r8_local(is:ie,js:je,npz+1)) ! U - if (is_master()) print*, 'Writing : ', TRIM(fname1), ' U' - r8_local(is:ie,js:je,1:npz) = FV_Atm(1)%u(is:ie,js:je,1:npz) - if (isBinFV) then - if (n_writers==1) then - call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),rc=status) - VERIFY_(status) - else - call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) - VERIFY_(status) - end if + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' U' + r8_local(is:ie,js:je,1:npz) = FV_Atm(1)%u(is:ie,js:je,1:npz) + if (isBinFV) then + if (n_writers==1) then + call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),rc=status) + VERIFY_(status) else - call MAPL_VarWrite(ncioOut,"U",r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) VERIFY_(status) - end if + end if + else + write(*,*)'bma writing u' + call MAPL_VarWrite(OutFmt,"U",r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + VERIFY_(status) + end if ! V - if (is_master()) print*, 'Writing : ', TRIM(fname1), ' V' - r8_local(is:ie,js:je,1:npz) = FV_Atm(1)%v(is:ie,js:je,1:npz) - if (isBinFV) then - if (n_writers==1) then - call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),rc=status) - VERIFY_(status) - else - call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) - VERIFY_(status) - end if + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' V' + r8_local(is:ie,js:je,1:npz) = FV_Atm(1)%v(is:ie,js:je,1:npz) + if (isBinFV) then + if (n_writers==1) then + call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),rc=status) + VERIFY_(status) else - call MAPL_VarWrite(ncioOut,"V",r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) - VERIFY_(status) - end if + call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + VERIFY_(status) + end if + else + call MAPL_VarWrite(OutFmt,"V",r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + VERIFY_(status) + end if ! PT - if (is_master()) print*, 'Writing : ', TRIM(fname1), ' PT' + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' PT' - if (isBinFV) then - if (n_writers==1) then - call MAPL_VarWrite(OUNIT,grid,pt_local(is:ie,js:je,1:npz),rc=status) - VERIFY_(status) - else - call MAPL_VarWrite(OUNIT,grid,pt_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) - VERIFY_(status) - end if + if (isBinFV) then + if (n_writers==1) then + call MAPL_VarWrite(OUNIT,grid,pt_local(is:ie,js:je,1:npz),rc=status) + VERIFY_(status) else - call MAPL_VarWrite(ncioOut,"PT",pt_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + call MAPL_VarWrite(OUNIT,grid,pt_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) VERIFY_(status) end if + else + call MAPL_VarWrite(OutFmt,"PT",pt_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + VERIFY_(status) + end if ! PE - if (is_master()) print*, 'Writing : ', TRIM(fname1), ' PE' - do k=1,npz+1 - r8_local(is:ie,js:je,k) = FV_Atm(1)%pe(is:ie,k,js:je) - enddo + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' PE' + do k=1,npz+1 + r8_local(is:ie,js:je,k) = FV_Atm(1)%pe(is:ie,k,js:je) + enddo + if (isBinFV) then + if (n_writers==1) then + call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz+1),rc=status) + VERIFY_(status) + else + call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz+1),arrdes=arrdes,rc=status) + VERIFY_(status) + end if + else + call MAPL_VarWrite(OutFmt,"PE",r8_local(is:ie,js:je,1:npz+1),arrdes=arrdes,rc=status) + VERIFY_(status) + end if +! PKZ + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' PKZ' + r8_local(is:ie,js:je,1:npz) = FV_Atm(1)%pkz(is:ie,js:je,1:npz) + if (isBinFV) then + if (n_writers==1) then + call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),rc=status) + VERIFY_(status) + else + call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + VERIFY_(status) + end if + else + call MAPL_VarWrite(OutFmt,"PKZ",r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + VERIFY_(status) + end if + + if (.not. fv_atm(1)%flagstruct%hydrostatic) then +! DZ + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' DZ' + r8_local(is:ie,js:je,1:npz) = FV_Atm(1)%delz(is:ie,js:je,1:npz) if (isBinFV) then if (n_writers==1) then - call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz+1),rc=status) + call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),rc=status) VERIFY_(status) else - call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz+1),arrdes=arrdes,rc=status) + call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) VERIFY_(status) end if else - call MAPL_VarWrite(ncioOut,"PE",r8_local(is:ie,js:je,1:npz+1),arrdes=arrdes,rc=status) + call MAPL_VarWrite(OutFmt,"DZ",r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) VERIFY_(status) end if -! PKZ - if (is_master()) print*, 'Writing : ', TRIM(fname1), ' PKZ' - r8_local(is:ie,js:je,1:npz) = FV_Atm(1)%pkz(is:ie,js:je,1:npz) + +! W + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' W' + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' DZ' + r8_local(is:ie,js:je,1:npz) = FV_Atm(1)%w(is:ie,js:je,1:npz) if (isBinFV) then if (n_writers==1) then call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),rc=status) @@ -605,62 +694,29 @@ program interp_restarts VERIFY_(status) end if else - call MAPL_VarWrite(ncioOut,"PKZ",r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + call MAPL_VarWrite(OutFmt,"W",r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) VERIFY_(status) end if + endif - if (.not. fv_atm(1)%flagstruct%hydrostatic) then -! DZ - if (is_master()) print*, 'Writing : ', TRIM(fname1), ' DZ' - r8_local(is:ie,js:je,1:npz) = FV_Atm(1)%delz(is:ie,js:je,1:npz) - if (isBinFV) then - if (n_writers==1) then - call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),rc=status) - VERIFY_(status) - else - call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) - VERIFY_(status) - end if - else - call MAPL_VarWrite(ncioOut,"DZ",r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) - VERIFY_(status) - end if - -! W - if (is_master()) print*, 'Writing : ', TRIM(fname1), ' W' - if (is_master()) print*, 'Writing : ', TRIM(fname1), ' DZ' - r8_local(is:ie,js:je,1:npz) = FV_Atm(1)%w(is:ie,js:je,1:npz) - if (isBinFV) then - if (n_writers==1) then - call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),rc=status) - VERIFY_(status) - else - call MAPL_VarWrite(OUNIT,grid,r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) - VERIFY_(status) - end if - else - call MAPL_VarWrite(ncioOut,"W",r8_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + if (isBinFV) then + if (n_writers > 1) then + if (AmWriter) then + call MPI_FILE_CLOSE(OUNIT,status) VERIFY_(status) end if - endif - - if (isBinFV) then - if (n_writers > 1) then - if (AmWriter) then - call MPI_FILE_CLOSE(OUNIT,status) - VERIFY_(status) - end if - else - close (OUNIT) - end if else - if (AmWriter) call MAPL_NCIOClose(ncioOut,destroy=.true.,rc=status) + close (OUNIT) end if + else + if (AmWriter) call OutFmt%close() + deallocate(InCfg,OutCfg) + end if - deallocate (r8_local) - deallocate (pt_local) + deallocate (r8_local) + deallocate (pt_local) - endif + endif ! MOIST ! @@ -681,73 +737,77 @@ program interp_restarts else ounit = getfile(trim(fname1),form='unformatted',rc=status) VERIFY_(status) - !open(OUNIT,file=fname1,access='sequential',form='unformatted') end if else - !if (AmWriter) then - ncio = MAPL_NCIOOpen("moist_internal_restart_in",rc=status) - imc = npx-1 - jmc = imc*6 - call MAPL_NCIOChangeRes(ncio,ncioOut,lonSize=imc,latSize=jmc,levSize=npz) - call MAPL_NCIOSet(ncioOut,filename=fname1) + imc = npx-1 + jmc = imc*6 + call InFmt%open("moist_internal_restart_in",pFIO_READ,rc=status) + allocate(InCfg(1),OutCfg(1)) + InCfg(1)=InFmt%read(rc=status) + call MAPL_IOChangeRes(InCfg(1),OutCfg(1),(/'lon','lat','lev'/),(/imc,jmc,npz/),rc=status) if (AmWriter) then - call MAPL_NCIOCreateFile(ncioOut,comm=arrdes%writers_comm,info=info,rc=status) + call OutFmt%create_par(fname1,comm=arrdes%writers_comm,info=info,rc=status) + call OutFmt%write(OutCfg(1),rc=status) + end if + deallocate(InCfg) + call InFmt%close() + end if + end if + ! binary path + if (isBinMoist) then + do iq=1,FV_Atm(1)%ncnst + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' ', iq + r4_local(is:ie,js:je,1:npz) = FV_Atm(1)%q(is:ie,js:je,:,iq) + iq0=iq + r4_local(is:ie,js:je,1:npz) = FV_Atm(1)%q(is:ie,js:je,:,iq0) + if (n_writers == 1) then + call MAPL_VarWrite(OUNIT,grid,r4_local(is:ie,js:je,1:npz),rc=status) VERIFY_(status) - call MAPL_NCIOClose(ncio,destroy=.true.,rc=status) + else + call MAPL_VarWrite(OUNIT,grid,r4_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + VERIFY_(status) + end if + end do + if (n_writers == 1) then + close (OUNIT) + else + if (AmWriter) then + call MPI_FILE_CLOSE(OUNIT,status) VERIFY_(status) end if end if + else + ! netcdf path + siter = all_moist_Vars%begin() + Variables => OutCfg(1)%get_variables() lcnt_var=2 - ! binary path - if (isBinMoist) then - do iq=1,FV_Atm(1)%ncnst - if (is_master()) print*, 'Writing : ', TRIM(fname1), ' ', iq - r4_local(is:ie,js:je,1:npz) = FV_Atm(1)%q(is:ie,js:je,:,iq) - iq0=iq - r4_local(is:ie,js:je,1:npz) = FV_Atm(1)%q(is:ie,js:je,:,iq0) - if (n_writers == 1) then - call MAPL_VarWrite(OUNIT,grid,r4_local(is:ie,js:je,1:npz),rc=status) - VERIFY_(status) + ivar=0 + do while (siter /= all_moist_vars%end()) + ivar=ivar+1 + var_name => siter%get() + myVariable => variables%at(var_name) + var_dimensions => myVariable%get_dimensions() + ndims = var_dimensions%size() + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' ', ivar + if (ndims==2) then + r4_local2d(is:ie,js:je)=0.0 + call MAPL_VarWrite(OutFmt,trim(var_name),r4_local2d(is:ie,js:je),arrdes=arrdes,rc=status) + VERIFY_(status) + else if (ndims==3) then + if (trim(var_name)=='Q') then + iq0=1 else - call MAPL_VarWrite(OUNIT,grid,r4_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) - VERIFY_(status) - end if - end do - if (n_writers == 1) then - close (OUNIT) - else - if (AmWriter) then - call MPI_FILE_CLOSE(OUNIT,status) - VERIFY_(status) + iq0=lcnt_var + lcnt_var=lcnt_var+1 end if + r4_local(is:ie,js:je,1:npz) = FV_Atm(1)%q(is:ie,js:je,:,iq0) + call MAPL_VarWrite(OutFmt,triM(var_name),r4_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) + VERIFY_(status) end if - else - ! netcdf path - call MAPL_NCIOGetDimSizes(ncioOut,nVars=nVars) - lcnt_var=2 - do ivar=1,nVars - call MAPL_NCIOGetVarName(ncioOut,ivar,vname) - if (is_master()) print*, 'Writing : ', TRIM(fname1), ' ', ivar - call MAPL_NCIOVarGetDims(ncioOut,vname,ndims,dimSize) - if (ndims==2) then - r4_local2d(is:ie,js:je)=0.0 - call MAPL_VarWrite(ncioOut,vname,r4_local2d(is:ie,js:je),arrdes=arrdes,rc=status) - VERIFY_(status) - else if (ndims==3) then - if (trim(vname)=='Q') then - iq0=1 - else - iq0=lcnt_var - lcnt_var=lcnt_var+1 - end if - r4_local(is:ie,js:je,1:npz) = FV_Atm(1)%q(is:ie,js:je,:,iq0) - call MAPL_VarWrite(ncioOut,vname,r4_local(is:ie,js:je,1:npz),arrdes=arrdes,rc=status) - VERIFY_(status) - end if - end do - if (AmWriter) call MAPL_NCIOClose(ncioOut,destroy=.true.,rc=status) - !if (AmWriter) call MAPL_NCIOClose(ncioOut,destroy=.true.,rc=status) - end if + call siter%next() + end do + if (AmWriter) call OutFmt%close() + deallocate(outCfg) end if deallocate(r4_local) @@ -774,20 +834,22 @@ program interp_restarts VERIFY_(status) end if else - !if (AmWriter) then - ncio = MAPL_NCIOOpen(trim(rst_files(ifile)%file_name),rc=status) + if (AmWriter) then imc = npx-1 jmc = imc*6 - call MAPL_NCIOChangeRes(ncio,ncioOut,lonSize=imc,latSize=jmc,levSize=npz) - call MAPL_NCIOSet(ncioOut,filename=fname1) - if (AmWriter) then - call MAPL_NCIOCreateFile(ncioOut,comm=arrdes%writers_comm,info=info,rc=status) + call InFmt%open(trim(rst_files(ifile)%file_name),pFIO_READ,rc=status) + allocate(InCfg(1),OutCfg(1)) + InCfg(1)=InFmt%read(rc=status) + call MAPL_IOChangeRes(InCfg(1),OutCfg(1),(/'lon','lat','lev'/),(/imc,jmc,npz/),rc=status) + call OutFmt%create_par(fname1,comm=arrdes%writers_comm,info=info,rc=status) + call OutFmt%write(OutCfg(1),rc=status) + deallocate(InCfg,OutCfg) + call InFmt%close() end if end if do iq=1,size(rst_files(ifile)%vars) - if (.not.rst_files(ifile)%isBin .and. AmWriter) then - call MAPL_NCIOGetVarName(ncio,iq,vname) - end if + vname = trim(rst_files(ifile)%vars(iq)%name) + if (is_master()) print*, 'Writing : ', TRIM(fname1), ' ', iq nlev = rst_files(ifile)%vars(iq)%nlev allocate(r4_local(is:ie,js:je,nlev)) if (rst_files(ifile)%isBin) then @@ -812,13 +874,11 @@ program interp_restarts end if else if (nlev/=1) then - if (is_master()) print*, 'Writing 3D: ', iq, TRIM(vname) r4_local(is:ie,js:je,1:nlev)=rst_files(ifile)%vars(iq)%ptr3d(is:ie,js:je,1:nlev) - call MAPL_VarWrite(ncioOut,vname,r4_local(is:ie,js:je,1:nlev),arrdes) + call MAPL_VarWrite(OutFmt,vname,r4_local(is:ie,js:je,1:nlev),arrdes) else - if (is_master()) print*, 'Writing 2D: ', TRIM(vname) r4_local2d(is:ie,js:je)=rst_files(ifile)%vars(iq)%ptr2d(is:ie,js:je) - call MAPL_VarWrite(ncioOut,vname,r4_local2d(is:ie,js:je),arrdes=arrdes) + call MAPL_VarWrite(OutFmt,vname,r4_local2d(is:ie,js:je),arrdes=arrdes) end if end if end do @@ -833,10 +893,8 @@ program interp_restarts end if end if else - !if (AmWriter) then - call MAPL_NCIOClose(ncio,destroy=.true.,rc=status) if (AmWriter) then - call MAPL_NCIOClose(ncioOut,destroy=.true.,rc=status) + call OutFmt%close() end if end if diff --git a/rs_scale.F90 b/rs_scale.F90 index 99ccaef..c5471cc 100644 --- a/rs_scale.F90 +++ b/rs_scale.F90 @@ -1,6 +1,5 @@ program main - use MAPL_IOMod - use MAPL_ConstantsMod, only: MAPL_PSDRY + use MAPL implicit none ! ************************************************************************ @@ -62,8 +61,8 @@ program main real*8, parameter :: eps = epsilon(1.0d-10) real :: kappa - type(MAPL_NCIO) :: InNCIODyn, OutNCIODyn - type(MAPL_NCIO) :: InNCIOMoist, OutNCIOMoist + type(Netcdf4_FileFormatter) :: InDyn,OutDyn,InMoist,OutMoist + type(FileMetadata) :: InCfgDyn,OutCfgDyn,InCfgMoist,OutCfgMoist integer :: filetypeDyn, filetypeMoist, nVars,nVarsMoist kappa = 2.0/7.0 @@ -92,11 +91,15 @@ program main if (filetypeDyn == 0) then - InNCIODyn = MAPL_NCIOOpen(dynrst,rc=rc) - InNCIOMoist = MAPL_NCIOOpen(mstrst,rc=rc) - call MAPL_NCIOGetDimSizes(InNCIODyn,lon=im,lat=jm,lev=lm,nVars=nVars) - call MAPL_NCIOGetDimSizes(InNCIOMoist,nVars=nVarsMoist) - + call InDyn%open(dynrst,pFIO_READ,rc=rc) + InCfgDyn=InDyn%read(rc=rc) + im = InCfgDyn%get_dimension('lon',rc=rc) + jm = InCfgDyn%get_dimension('lat',rc=rc) + lm = InCfgDyn%get_dimension('lev',rc=rc) + call MAPL_IOCountNonDimVars(InCfgDyn,nVars) + call InMoist%open(mstrst,pFIO_READ,rc=rc) + InCfgMoist = InMoist%read(rc=rc) + call MAPL_IOCountNonDimVars(InCfgMoist,nVarsMoist) nymd=0 nhms=0 @@ -144,32 +147,32 @@ program main if (filetypeDyn == 0) then - call MAPL_VarRead(InNCIODyn,"AK",ak) - call MAPL_VarRead(InNCIODyn,"BK",bk) + call MAPL_VarRead(InDyn,"AK",ak) + call MAPL_VarRead(InDyn,"BK",bk) do L=1,lm - call MAPL_VarRead(InNCIODyn,"U",u(:,:,L),lev=l) + call MAPL_VarRead(InDyn,"U",u(:,:,L),lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIODyn,"V",v(:,:,L),lev=l) + call MAPL_VarRead(InDyn,"V",v(:,:,L),lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIODyn,"PT",th(:,:,L),lev=l) + call MAPL_VarRead(InDyn,"PT",th(:,:,L),lev=l) enddo do L=1,lm+1 - call MAPL_VarRead(InNCIODyn,"PE",ple(:,:,L),lev=l) + call MAPL_VarRead(InDyn,"PE",ple(:,:,L),lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIODyn,"PKZ",pk(:,:,L),lev=l) + call MAPL_VarRead(InDyn,"PKZ",pk(:,:,L),lev=l) enddo ! check if we have delz and w if (nvars == 9) then allocate( delz(im,jm,lm ) ) allocate( w(im,jm,lm ) ) do L=1,lm - call MAPL_VarRead(InNCIODyn,"DZ",delz(:,:,L),lev=l) + call MAPL_VarRead(InDyn,"DZ",delz(:,:,L),lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIODyn,"W",w(:,:,L),lev=l) + call MAPL_VarRead(InDyn,"W",w(:,:,L),lev=l) enddo end if @@ -218,25 +221,25 @@ program main if (filetypeMoist == 0) then do L=1,lm - call MAPL_VarRead(InNCIOMoist,"Q",qv(:,:,l),lev=l) + call MAPL_VarRead(InMoist,"Q",qv(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIOMoist,"QLLS",qlls(:,:,l),lev=l) + call MAPL_VarRead(InMoist,"QLLS",qlls(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIOMoist,"QLCN",qlcn(:,:,l),lev=l) + call MAPL_VarRead(InMoist,"QLCN",qlcn(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIOMoist,"CLLS",cfls(:,:,l),lev=l) + call MAPL_VarRead(InMoist,"CLLS",cfls(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIOMoist,"CLCN",cfcn(:,:,l),lev=l) + call MAPL_VarRead(InMoist,"CLCN",cfcn(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIOMoist,"QILS",qils(:,:,l),lev=l) + call MAPL_VarRead(InMoist,"QILS",qils(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIOMoist,"QICN",qicn(:,:,l),lev=l) + call MAPL_VarRead(InMoist,"QICN",qicn(:,:,l),lev=l) enddo else @@ -384,33 +387,33 @@ program main if (filetypeDyn == 0) then dynrst = trim(dynrst) // '.scaled' - call MAPL_NCIOChangeRes(InNCIODyn,OutNCIODyn,latSize=jm,lonSize=im) - call MAPL_NCIOSet(OutNCIODyn,filename=dynrst) - call MAPL_NCIOCreateFile(OutNCIODyn) + OutCfgDyn=InCfgDyn + call OutDyn%create(dynrst,rc=rc) + call OutDyn%write(OutCfgDyn,rc=rc) - call MAPL_VarWrite(OutNCIODyn,"AK",ak) - call MAPL_VarWrite(OutNCIODyn,"BK",bk) + call MAPL_VarWrite(OutDyn,"AK",ak) + call MAPL_VarWrite(OutDyn,"BK",bk) do L=1,lm - call MAPL_VarWrite(OutNCIODyn,"U",u(:,:,L),lev=L) + call MAPL_VarWrite(OutDyn,"U",u(:,:,L),lev=L) enddo do L=1,lm - call MAPL_VarWrite(OutNCIODyn,"V",v(:,:,L),lev=L) + call MAPL_VarWrite(OutDyn,"V",v(:,:,L),lev=L) enddo do L=1,lm - call MAPL_VarWrite(OutNCIODyn,"PT",th(:,:,L),lev=L) + call MAPL_VarWrite(OutDyn,"PT",th(:,:,L),lev=L) enddo do L=1,lm+1 - call MAPL_VarWrite(OutNCIODyn,"PE",ple(:,:,L),lev=L) + call MAPL_VarWrite(OutDyn,"PE",ple(:,:,L),lev=L) enddo do L=1,lm - call MAPL_VarWrite(OutNCIODyn,"PKZ",pk(:,:,L),lev=L) + call MAPL_VarWrite(OutDyn,"PKZ",pk(:,:,L),lev=L) enddo if (nvars==9) then do L=1,lm - call MAPL_VarWrite(OutNCIODyn,"DZ",delz(:,:,L),lev=L) + call MAPL_VarWrite(OutDyn,"DZ",delz(:,:,L),lev=L) enddo do L=1,lm - call MAPL_VarWrite(OutNCIODyn,"W",w(:,:,L),lev=L) + call MAPL_VarWrite(OutDyn,"W",w(:,:,L),lev=L) enddo end if @@ -474,42 +477,43 @@ program main if (fileTypeMoist == 0) then mstrst = trim(mstrst) // '.scaled' - call MAPL_NCIOChangeRes(InNCIOMoist,OutNCIOMoist,latSize=jm,lonSize=im) - call MAPL_NCIOSet(OutNCIOMoist,filename=mstrst) - call MAPL_NCIOCreateFile(OutNCIOMoist) + + OutCfgMoist=InCfgMoist + call OutMoist%create(mstrst,rc=rc) + call OutMoist%write(OutCfgMoist,rc=rc) do L=1,lm - call MAPL_VarWrite(OutNCIOMoist,"Q",qv(:,:,l),lev=l) + call MAPL_VarWrite(OutMoist,"Q",qv(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarWrite(OutNCIOMoist,"QLLS",qlls(:,:,l),lev=l) + call MAPL_VarWrite(OutMoist,"QLLS",qlls(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarWrite(OutNCIOMoist,"QLCN",qlcn(:,:,l),lev=l) + call MAPL_VarWrite(OutMoist,"QLCN",qlcn(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarWrite(OutNCIOMoist,"CLLS",cfls(:,:,l),lev=l) + call MAPL_VarWrite(OutMoist,"CLLS",cfls(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarWrite(OutNCIOMoist,"CLCN",cfcn(:,:,l),lev=l) + call MAPL_VarWrite(OutMoist,"CLCN",cfcn(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarWrite(OutNCIOMoist,"QILS",qils(:,:,l),lev=l) + call MAPL_VarWrite(OutMoist,"QILS",qils(:,:,l),lev=l) enddo do L=1,lm - call MAPL_VarWrite(OutNCIOMoist,"QICN",qicn(:,:,l),lev=l) + call MAPL_VarWrite(OutMoist,"QICN",qicn(:,:,l),lev=l) enddo ! check if we have ncpl and ncpi if (nVarsMoist == 9) then allocate( dum4(im,jm) ) do L=1,lm - call MAPL_VarRead(InNCIOMoist,"NCPL",dum4,lev=l) - call MAPL_VarWrite(OutNCIOMoist,"NCPL",dum4,lev=l) + call MAPL_VarRead(InMoist,"NCPL",dum4,lev=l) + call MAPL_VarWrite(OutMoist,"NCPL",dum4,lev=l) enddo do L=1,lm - call MAPL_VarRead(InNCIOMoist,"NCPI",dum4,lev=l) - call MAPL_VarWrite(OutNCIOMoist,"NCPI",dum4,lev=l) + call MAPL_VarRead(InMoist,"NCPI",dum4,lev=l) + call MAPL_VarWrite(OutMoist,"NCPI",dum4,lev=l) enddo end if