diff --git a/CMakeLists.txt b/CMakeLists.txt index 282da6c3b..57b335525 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -56,7 +56,7 @@ endif () esma_add_library (${this} SRCS ${srcs} - DEPENDENCIES MAPL_Base + DEPENDENCIES MAPL gftl-shared INCLUDES ${extra_incs} ${INC_ESMF} ${INC_NETCDF}) if (FV_PRECISION STREQUAL R4) diff --git a/geos_utils/fv_regrid_c2c.F90 b/geos_utils/fv_regrid_c2c.F90 index c7437ca85..ad5df1eee 100644 --- a/geos_utils/fv_regrid_c2c.F90 +++ b/geos_utils/fv_regrid_c2c.F90 @@ -1,7 +1,7 @@ module fv_regrid_c2c #ifdef MAPL_MODE -#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 +#define DEALLOCGLOB_(A) if(associated(A)) then;A=0;if(MAPL_ShmInitialized) then; call MAPL_DeAllocNodeArray(A,rc=status);else; deallocate(A);endif;NULLIFY(A);endif #endif use fv_arrays_mod, only: REAL4, REAL8, FVPRC @@ -13,11 +13,9 @@ module fv_regrid_c2c use tracer_manager_mod, only: get_tracer_names, get_number_tracers, get_tracer_index use field_manager_mod, only: MODEL_ATMOS - use MAPL_MOD, only: MAPL_PI_R8, MAPL_OMEGA, MAPL_GRAV, & - MAPL_KAPPA, MAPL_RGAS, MAPL_RVAP, & - MAPL_CP - use MAPL_IOMod - use MAPL_ShmemMod + use MAPL + use gFTL_StringVector + use gFTL_StringIntegerMap use, intrinsic :: iso_fortran_env, only: REAL64, REAL32 use fv_arrays_mod, only: fv_atmos_type, fv_grid_type, fv_grid_bounds_type, FVPRC, REAL4, REAL8 @@ -34,8 +32,7 @@ module fv_regrid_c2c use mpp_mod, only: mpp_pe use memutils_mod, only: print_memuse_stats use fv_regridding_utils - use pFIO_StringVectorMod - use pFIO_StringIntegerMapMod + use ESMF implicit none @@ -63,11 +60,12 @@ module fv_regrid_c2c contains - subroutine get_geos_ic( Atm, extra_rst, rstcube) + subroutine get_geos_ic( Atm, extra_rst, rstcube, gridOut) type(fv_atmos_type), intent(inout) :: Atm(:) type(fv_rst), pointer, intent(inout) :: extra_rst(:) logical :: rstcube + type(ESMF_Grid), intent(inout) :: gridOut real(FVPRC):: alpha = 0. integer i,j logical :: cubed_sphere,fv_diag_ic @@ -99,14 +97,14 @@ subroutine get_geos_ic( Atm, extra_rst, rstcube) if (extra_rst(i)%have_descriptor) then do j=1,size(extra_rst(i)%vars) if (extra_rst(i)%vars(j)%nLev/=1) then - allocate(extra_rst(i)%vars(j)%ptr3d(isd:ied,jsd:jed,extra_rst(i)%vars(j)%nLev),source=0.0d0 ) + allocate(extra_rst(i)%vars(j)%ptr3d(isd:ied,jsd:jed,extra_rst(i)%vars(j)%nLev),source=0.0 ) else - allocate(extra_rst(i)%vars(j)%ptr2d(isd:ied,jsd:jed), source=0.0d0 ) + allocate(extra_rst(i)%vars(j)%ptr2d(isd:ied,jsd:jed), source=0.0 ) end if enddo else do j=1,size(extra_rst(i)%vars) - allocate(extra_rst(i)%vars(j)%ptr3d(isd:ied,jsd:jed,extra_rst(i)%vars(j)%nLev),source=0.0d0 ) + allocate(extra_rst(i)%vars(j)%ptr3d(isd:ied,jsd:jed,extra_rst(i)%vars(j)%nLev),source=0.0 ) enddo end if enddo @@ -118,7 +116,7 @@ subroutine get_geos_ic( Atm, extra_rst, rstcube) else if (allocated(Atm(1)%q)) deallocate( Atm(1)%q ) allocate ( Atm(1)%q(isd:ied,jsd:jed,Atm(1)%npz,Atm(1)%ncnst) ) - call get_geos_cubed_ic( Atm, extra_rst ) + call get_geos_cubed_ic( Atm, extra_rst, gridOut ) endif call prt_maxmin('T', Atm(1)%pt, is, ie, js, je, ng, Atm(1)%npz, 1.0_FVPRC) @@ -132,12 +130,13 @@ subroutine get_geos_ic( Atm, extra_rst, rstcube) end subroutine get_geos_ic - subroutine get_geos_cubed_ic( Atm, extra_rst ) + subroutine get_geos_cubed_ic( Atm, extra_rst, gridOut ) use GHOST_CUBSPH_mod, only : A_grid, ghost_cubsph_update use CUB2CUB_mod, only: get_c2c_weight, & interpolate_c2c type(fv_atmos_type), intent(inout) :: Atm(:) type(fv_rst), pointer, intent(inout) :: extra_rst(:) + type(ESMF_Grid), intent(inout) :: gridOut character(len=128) :: fname, fname1 real(FVPRC), allocatable:: pkz0(:,:) @@ -157,8 +156,7 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) real(REAL64), allocatable :: akbk_r8(:) - real(REAL64), dimension(:,:,:,:), allocatable :: corner_in, corner_out, weight_c2c - integer, dimension(:,:,:,:), allocatable :: index_c2c + real(REAL64), dimension(:,:,:,:), allocatable :: corner_in, corner_out integer :: is_i,ie_i, js_i,je_i integer :: isd_i,ied_i, jsd_i,jed_i @@ -175,21 +173,31 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) integer :: filetype,nqmap logical :: isNC4 - type(MAPL_NCIO) :: ncio + type(Netcdf4_Fileformatter) :: formatter + type(FileMetadata), allocatable :: cfg(:) integer :: nDims, nVars, ivar, dimSizes(3) character(len=128) :: vname real(FVPRC), allocatable :: gslice_r4(:,:) real*8, allocatable :: gslice_r8(:,:) - integer :: tileoff,lvar_cnt,ifile,nlev,dimSize(3) + integer :: tileoff,lvar_cnt,ifile,nlev type(fv_rst), pointer :: tracer_bundles(:) => null() !bma added + character(len=:), pointer :: var_name + type(StringVariableMap), pointer :: variables + type(Variable), pointer :: myVariable + type(StringVector) :: all_moist_vars + type(StringVector), pointer :: var_dimensions + type(StringVectorIterator) :: siter type(StringVector) :: moist_variables type(StringIntegerMap) :: moist_tracers integer, pointer :: iptr character(len=:), pointer :: cptr type(StringIntegerMapIterator) :: iter character(len=128) :: moist_order(9) = (/"Q ","QLLS","QLCN","CLLS","CLCN","QILS","QICN","NCPL","NCPI"/) + type(CubedSphereGridFactory) :: cs_factory + type(ESMF_Grid) :: gridIn + class(AbstractRegridder), pointer :: regridder=>null() npx = Atm(1)%npx npy = Atm(1)%npy @@ -212,9 +220,13 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) if (isNC4) then - NCIO = MAPL_NCIOOpen(fname) - call MAPL_NCIOGetDimSizes(NCIO,lon=im,lat=jm,lev=km) - allocate(gslice_r8(im,jm)) + allocate(cfg(1)) + call formatter%open(fname,pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + im =cfg(1)%get_dimension('lon',rc=status) + jm =cfg(1)%get_dimension('lat',rc=status) + km =cfg(1)%get_dimension('lev',rc=status) + allocate(gslice_r8(im,jm),stat=status) else @@ -240,6 +252,9 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) & field not found') endif + cs_factory = CubedSphereGridFactory(nx=Atm(1)%layout(1),ny=Atm(1)%layout(2),im_world=im,lm=km) + gridIn = grid_manager%make_grid(cs_factory,rc=status) + regridder => new_regridder_manager%make_regridder(gridIn,gridOut,REGRID_METHOD_BILINEAR,rc=status) !--------------------------------------------------------------------! ! setup input cubed-sphere domain ! !--------------------------------------------------------------------! @@ -264,18 +279,6 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) enddo enddo call print_memuse_stats('get_geos_cubed_ic: init corner_out') -!--------------------------------------------------------------------! -! calculate weights and indices from bilinear interpolation ! -! from grid_in to grid_out ! -!--------------------------------------------------------------------! - allocate(index_c2c (3, is:ie, js:je, tile:tile)) - allocate(weight_c2c(4, is:ie, js:je, tile:tile)) - call get_c2c_weight(ntiles, im+1, (jm/6)+1, & - is_i, ie_i, js_i, je_i, isd_i, ied_i, jsd_i, jed_i, & - corner_in(:,is_i:ie_i+1,js_i:je_i+1,tile), & - npx, npy, is,ie, js,je, tile,tile, & - corner_out(:,is:ie+1,js:je+1,tile:tile), & - index_c2c, weight_c2c, domain_i) npts = im call print_memuse_stats('get_geos_cubed_ic: get_c2c_weight') @@ -283,13 +286,13 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) allocate ( bk0(km+1) ) allocate ( akbk_r8(km+1) ) if (isNC4) then - call MAPL_VarRead(NCIO,"AK",akbk_r8) + call MAPL_VarRead(formatter,"AK",akbk_r8) else read (IUNIT, IOSTAT=status) akbk_r8 end if ak0 = akbk_r8 if (isNC4) then - call MAPL_VarRead(NCIO,"BK",akbk_r8) + call MAPL_VarRead(formatter,"BK",akbk_r8) else read (IUNIT, IOSTAT=status) akbk_r8 end if @@ -304,7 +307,7 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) if (isNC4) then tileoff = (tile-1)*(jm/ntiles) do k=1,km - call MAPL_VarRead(NCIO,"U",gslice_r8,lev=k) + call MAPL_VarRead(formatter,"U",gslice_r8,lev=k) u0(is_i:ie_i,js_i:je_i,k) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i) enddo else @@ -320,7 +323,7 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) if (isNC4) then tileoff = (tile-1)*(jm/ntiles) do k=1,km - call MAPL_VarRead(NCIO,"V",gslice_r8,lev=k) + call MAPL_VarRead(formatter,"V",gslice_r8,lev=k) v0(is_i:ie_i,js_i:je_i,k) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i) enddo else @@ -355,7 +358,7 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) allocate ( va(is:ie,js:je,km) ) call interp_c2c_vect(npts, npts, npx-1, npy-1, km, ntiles, domain_i, & is,ie, js,je, isd_i,ied_i, jsd_i,jed_i, is_i,ie_i, js_i,je_i, & - u0, v0, ua, va, index_c2c, weight_c2c, corner_in(:,is_i:ie_i+1,js_i:je_i+1,tile), corner_out, Atm(1)%gridstruct) + u0, v0, ua, va, regridder, corner_in(:,is_i:ie_i+1,js_i:je_i+1,tile), corner_out, Atm(1)%gridstruct) deallocate ( v0 ) deallocate ( u0 ) deallocate ( corner_in ) @@ -366,7 +369,7 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) if (isNC4) then tileoff = (tile-1)*(jm/ntiles) do k=1,km - call MAPL_VarRead(NCIO,"PT",gslice_r8,lev=k) + call MAPL_VarRead(formatter,"PT",gslice_r8,lev=k) t0(is_i:ie_i,js_i:je_i,k) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i) enddo else @@ -379,7 +382,7 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) ps0(:,:) = 0.0 if (isNC4) then tileoff = (tile-1)*(jm/ntiles) - call MAPL_VarRead(NCIO,"PE",gslice_r8,lev=km+1) + call MAPL_VarRead(formatter,"PE",gslice_r8,lev=km+1) ps0(is_i:ie_i,js_i:je_i) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i) else slice_2d = npts*npts*ntiles @@ -396,7 +399,7 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) if (isNC4) then tileoff = (tile-1)*(jm/ntiles) do k=1,km - call MAPL_VarRead(NCIO,"PKZ",gslice_r8,lev=k) + call MAPL_VarRead(formatter,"PKZ",gslice_r8,lev=k) pkz0(is_i:ie_i,js_i:je_i) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i) t0(is_i:ie_i,js_i:je_i,k) = t0(is_i:ie_i,js_i:je_i,k)*pkz0(is_i:ie_i,js_i:je_i) enddo @@ -412,7 +415,8 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) deallocate ( pkz0 ) if (isNC4) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) deallocate(gslice_r8) end if @@ -443,29 +447,11 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) ! Horiz Interp for surface pressure call prt_maxmin('PS_geos', ps0, is_i, ie_i, js_i, je_i, ng_i, 1, 1.0_FVPRC) - do j=js,je - do i=is,ie - ic=index_c2c(1,i,j,tile) - jc=index_c2c(2,i,j,tile) - psc(i,j)=weight_c2c(1,i,j,tile)*ps0(ic ,jc ) & - +weight_c2c(2,i,j,tile)*ps0(ic ,jc+1) & - +weight_c2c(3,i,j,tile)*ps0(ic+1,jc+1) & - +weight_c2c(4,i,j,tile)*ps0(ic+1,jc ) - enddo - enddo + call regridder%regrid(ps0(is_i:ie_i,js_i:je_i),psc(is:ie,js:je),rc=status) deallocate ( ps0 ) ! Horiz Interp for surface height call prt_maxmin('GZ_geos', gz0, is_i, ie_i, js_i, je_i, ng_i, 1, 1.0/grav) - do j=js,je - do i=is,ie - ic=index_c2c(1,i,j,tile) - jc=index_c2c(2,i,j,tile) - gzc(i,j)=weight_c2c(1,i,j,tile)*gz0(ic ,jc ) & - +weight_c2c(2,i,j,tile)*gz0(ic ,jc+1) & - +weight_c2c(3,i,j,tile)*gz0(ic+1,jc+1) & - +weight_c2c(4,i,j,tile)*gz0(ic+1,jc ) - enddo - enddo + call regridder%regrid(gz0(is_i:ie_i,js_i:je_i),gzc(is:ie,js:je),rc=status) deallocate ( gz0 ) ! Horiz Interp for Q @@ -485,21 +471,29 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) if (filetype /= 0) then offset=4 else + lvar_cnt = 0 allocate(gslice_r4(im,jm)) - NCIO = MAPL_NCIOOpen("moist_internal_restart_in") - call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars) - do ivar=1,nVars - call MAPL_NCIOGetVarName(NCIO,ivar,vname) - call MAPL_NCIOVarGetDims(NCIO,vname,ndims,dimSize) - if (ndims==3) call moist_variables%push_back(trim(vname)) !only do the 3d-variables + allocate(cfg(1)) + call formatter%open("moist_internal_restart_in",pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + all_moist_vars = MAPL_IOGetNonDimVars(cfg(1)) + siter = all_moist_vars%begin() + Variables => cfg(1)%get_variables() + 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 == 3) call moist_variables%push_back(trim(var_name)) + call siter%next() enddo - lvar_cnt=2 if (moist_variables%size() /= atm(1)%ncnst) call mpp_error(FATAL,'Wrong number of variables in moist file') tileoff = (tile-1)*(jm/ntiles) + lvar_cnt=2 end if do ivar=1,Atm(1)%ncnst - if (filetype /=0) then + if (filetype /= 0) then call moist_tracers%insert(trim(moist_order(ivar)),ivar) iq=ivar else @@ -517,27 +511,17 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) call parallel_read_file_r4('moist_internal_restart_in', npts, is_i,ie_i, js_i,je_i, 1, offset, q0(is_i:ie_i,js_i:je_i,k)) call mpp_update_domains(q0(:,:,k), domain_i) else - call MAPL_VarRead(NCIO,vname,gslice_r4,lev=k) + call MAPL_VarRead(formatter,vname,gslice_r4,lev=k) q0(is_i:ie_i,js_i:je_i,k)=gslice_r4(is_i:ie+i,tileoff+js_i:tileoff+je_i) end if - call mpp_update_domains(q0(:,:,k), domain_i) - do j=js,je - do i=is,ie - ic=index_c2c(1,i,j,tile) - jc=index_c2c(2,i,j,tile) - qp(i,j,k,iq)=weight_c2c(1,i,j,tile)*q0(ic ,jc ,k) & - +weight_c2c(2,i,j,tile)*q0(ic ,jc+1,k) & - +weight_c2c(3,i,j,tile)*q0(ic+1,jc+1,k) & - +weight_c2c(4,i,j,tile)*q0(ic+1,jc ,k) - enddo - enddo - + call regridder%regrid(q0(is_i:ie_i,js_i:je_i,k),qp(:,:,k,iq),rc=status) enddo call prt_maxmin( 'Q_geos_moist', q0, is_i, ie_i, js_i, je_i, ng_i, km, 1._FVPRC) enddo if (filetype == 0) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) deallocate(gslice_r4) end if @@ -576,8 +560,10 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) offset=4 else allocate(gslice_r4(im,jm)) - NCIO = MAPL_NCIOOpen(trim(tracer_bundles(ifile)%file_name)) - call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars) + allocate(cfg(1)) + call formatter%open(trim(tracer_bundles(ifile)%file_name),pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + call MAPL_IOCountNonDimVars(cfg(1),nvars,rc=status) tileoff = (tile-1)*(jm/ntiles) end if @@ -586,47 +572,30 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) do ivar=1,size(tracer_bundles(ifile)%vars) nlev=tracer_bundles(ifile)%vars(ivar)%nLev - if (filetype == 0) call MAPL_NCIOGetVarName(NCIO,ivar,vname) + if (filetype == 0) vname = trim(tracer_bundles(ifile)%vars(ivar)%name) do k=1,nlev if (filetype /= 0) then call parallel_read_file_r4(trim(tracer_bundles(ifile)%file_name), npts, is_i,ie_i, js_i,je_i, 1, offset, qlev(is_i:ie_i,js_i:je_i)) else if (tracer_bundles(ifile)%vars(ivar)%nLev/=1) then - call MAPL_VarRead(NCIO,vname,gslice_r4,lev=k) + call MAPL_VarRead(formatter,vname,gslice_r4,lev=k) else - call MAPL_VarRead(NCIO,vname,gslice_r4) + call MAPL_VarRead(formatter,vname,gslice_r4) end if qlev(is_i:ie_i,js_i:je_i)=gslice_r4(is_i:ie+i,tileoff+js_i:tileoff+je_i) end if - call mpp_update_domains(qlev, domain_i) - do j=js,je - do i=is,ie - ic=index_c2c(1,i,j,tile) - jc=index_c2c(2,i,j,tile) - if (tracer_bundles(ifile)%vars(ivar)%nLev/=1) then - tracer_bundles(ifile)%vars(ivar)%ptr3d(i,j,k) & - =weight_c2c(1,i,j,tile)*qlev(ic ,jc) & - +weight_c2c(2,i,j,tile)*qlev(ic ,jc+1) & - +weight_c2c(3,i,j,tile)*qlev(ic+1,jc+1) & - +weight_c2c(4,i,j,tile)*qlev(ic+1,jc) - else - tracer_bundles(ifile)%vars(ivar)%ptr2d(i,j) & - =weight_c2c(1,i,j,tile)*qlev(ic ,jc) & - +weight_c2c(2,i,j,tile)*qlev(ic ,jc+1) & - +weight_c2c(3,i,j,tile)*qlev(ic+1,jc+1) & - +weight_c2c(4,i,j,tile)*qlev(ic+1,jc) - end if - enddo - enddo - - - + if (tracer_bundles(ifile)%vars(ivar)%nLev/=1) then + call regridder%regrid(qlev(is_i:ie_i,js_i:js_i),tracer_bundles(ifile)%vars(ivar)%ptr3d(is:ie,js:je,k),rc=status) + else + call regridder%regrid(qlev(is_i:ie_i,js_i:js_i),tracer_bundles(ifile)%vars(ivar)%ptr2d(is:ie,js:je),rc=status) + end if enddo !call prt_maxmin( 'Q_geos_gocart', q0, is_i, ie_i, js_i, je_i, ng_i, km, 1._FVPRC) enddo if (filetype == 0) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) deallocate(gslice_r4) end if deallocate(qlev) @@ -639,20 +608,9 @@ subroutine get_geos_cubed_ic( Atm, extra_rst ) call prt_maxmin( 'T_geos', t0, is_i, ie_i, js_i, je_i, ng_i, km, 1.0_FVPRC) allocate ( tp(is:ie,js:je,km) ) do k=1,km - do j=js,je - do i=is,ie - ic=index_c2c(1,i,j,tile) - jc=index_c2c(2,i,j,tile) - tp(i,j,k)=weight_c2c(1,i,j,tile)*t0(ic ,jc ,k) & - +weight_c2c(2,i,j,tile)*t0(ic ,jc+1,k) & - +weight_c2c(3,i,j,tile)*t0(ic+1,jc+1,k) & - +weight_c2c(4,i,j,tile)*t0(ic+1,jc ,k) - enddo - enddo + call regridder%regrid(t0(is_i:ie_i,js_i:je_i,k),tp(:,:,k),rc=status) enddo deallocate ( t0 ) - deallocate( index_c2c ) - deallocate( weight_c2c ) ! Horz/Vert remap for scalars nqmap = Atm(1)%ncnst @@ -734,7 +692,6 @@ function match(var_name) result(inList) if (trim(exclude_vars(i))==trim(var_name)) inList = .true. enddo end function match - end subroutine get_geos_cubed_ic @@ -747,7 +704,7 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) real(FVPRC), allocatable:: ps0(:,:), gz0(:,:), t0(:,:,:), q0(:,:,:) real(FVPRC), allocatable:: u0(:,:,:), v0(:,:,:), ua0(:,:), va0(:,:) real(FVPRC), allocatable:: lat(:), lon(:), ak0(:), bk0(:) - integer :: i, j, k, l, iq, j1, j2, im, jm, km, npz + integer :: i, j, k, iq, j1, j2, im, jm, km, npz integer :: header(6) character (len=8) :: imc, jmc @@ -766,12 +723,13 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) integer :: filetype logical :: isNC4 - type(MAPL_NCIO) :: ncio + type(Netcdf4_Fileformatter) :: formatter + type(FileMetadata), allocatable :: cfg(:) integer :: nDims, nVars, ivar, dimSizes(3) character(len=128) :: vname integer :: lvar_cnt !bma added - type(stringVector) :: moist_variables + character(len=128) :: moist_order(9) = (/"Q ","QLLS","QLCN","CLLS","CLCN","QILS","QICN","NCPL","NCPI"/) type(fv_rst), pointer :: tracer_bundles(:) => null() integer :: ifile,nlev real(FVPRC), allocatable :: gslice_r4(:,:) @@ -796,8 +754,12 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) if (isNC4) then - NCIO = MAPL_NCIOOpen(fname) - call MAPL_NCIOGetDimSizes(NCIO,lon=im,lat=jm,lev=km) + allocate(cfg(1)) + call formatter%open(fname,pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + im =cfg(1)%get_dimension('lon',rc=status) + jm =cfg(1)%get_dimension('lat',rc=status) + km =cfg(1)%get_dimension('lev',rc=status) else @@ -831,13 +793,13 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) allocate ( bk0(km+1) ) allocate ( akbk_r8(km+1) ) if (isNC4) then - call MAPL_VarRead(NCIO,"AK",akbk_r8) + call MAPL_VarRead(formatter,"AK",akbk_r8) else read (IUNIT, IOSTAT=status) akbk_r8 end if ak0 = akbk_r8 if (isNC4) then - call MAPL_VarRead(NCIO,"BK",akbk_r8) + call MAPL_VarRead(formatter,"BK",akbk_r8) else read (IUNIT, IOSTAT=status) akbk_r8 end if @@ -851,7 +813,7 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) allocate ( u0(im,jm,km) ) do k=1,km if (isNC4) then - call MAPL_VarRead(NCIO,"U",r8latlon,lev=k) + call MAPL_VarRead(formatter,"U",r8latlon,lev=k) else read (IUNIT, IOSTAT=status) r8latlon end if @@ -864,7 +826,7 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) allocate ( v0(im,jm,km) ) do k=1,km if (isNC4) then - call MAPL_VarRead(NCIO,"V",r8latlon,lev=k) + call MAPL_VarRead(formatter,"V",r8latlon,lev=k) else read (IUNIT, IOSTAT=status) r8latlon end if @@ -912,7 +874,7 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) allocate ( t0(im,jm,km) ) do k=1,km if (isNC4) then - call MAPL_VarRead(NCIO,"PT",r8latlon,lev=k) + call MAPL_VarRead(formatter,"PT",r8latlon,lev=k) else read (IUNIT, IOSTAT=status) r8latlon end if @@ -924,7 +886,7 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) ! Read PE do k=1,km+1 if (isNC4) then - call MAPL_VarRead(NCIO,"PE",r8latlon,lev=k) + call MAPL_VarRead(formatter,"PE",r8latlon,lev=k) else read (IUNIT, IOSTAT=status) r8latlon end if @@ -936,7 +898,7 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) allocate ( pkz0(im,jm) ) do k=1,km if (isNC4) then - call MAPL_VarRead(NCIO,"PKZ",r8latlon,lev=k) + call MAPL_VarRead(formatter,"PKZ",r8latlon,lev=k) else read (IUNIT, IOSTAT=status) r8latlon end if @@ -950,7 +912,8 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) call print_memuse_stats('get_geos_latlon_ic: converted T') deallocate ( pkz0 ) if (isNC4) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) else close (IUNIT) end if @@ -1039,33 +1002,20 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) open(IUNIT,file="moist_internal_restart_in" ,access='sequential',form='unformatted',status='old') else lvar_cnt = 0 - NCIO = MAPL_NCIOOpen("moist_internal_restart_in") - call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars) - do ivar=1,nVars - call MAPL_NCIOGetVarName(NCIO,ivar,vname) - call moist_variables%push_back(trim(vname)) - enddo - lvar_cnt=2 + call formatter%open("moist_internal_restart_in",pFIO_READ,rc=status) + cfg = formatter%read(rc=status) + call MAPL_IOCountNonDimVars(cfg(1),nvars,rc=status) if (nVars /= Atm(1)%ncnst) call mpp_error(FATAL,'Wrong number of variables in moist file') end if do ivar=1,Atm(1)%ncnst - if (filetype /=0) then - iq=ivar - else - vname = moist_variables%at(ivar) - if (trim(vname)=='Q') then - iq=1 - else - iq=lvar_cnt - lvar_cnt=lvar_cnt+1 - end if - end if + if (filetype ==0) lvar_cnt=lvar_cnt+1 do k=1,km if (filetype /= 0) then read (IUNIT, IOSTAT=status) r4latlon else - call MAPL_VarRead(NCIO,vname,r4latlon,lev=k) + vname = trim(moist_order(lvar_cnt)) + call MAPL_VarRead(formatter,vname,r4latlon,lev=k) end if q0(1 :im/2,:,k) = r4latlon(im/2 + 1 :im , :) ! Regrid from -180:180 to 0:360 q0(im/2 + 1:im ,:,k) = r4latlon(1 :im/2, :) ! Regrid from -180:180 to 0:360 @@ -1074,7 +1024,7 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) i1 = id1(i,j) i2 = id2(i,j) j1 = jdc(i,j) - qp(i,j,k,iq) = s2c(i,j,1)*q0(i1,j1 ,k) + s2c(i,j,2)*q0(i2,j1 ,k) + & + qp(i,j,k,ivar) = s2c(i,j,1)*q0(i1,j1 ,k) + s2c(i,j,2)*q0(i2,j1 ,k) + & s2c(i,j,3)*q0(i2,j1+1,k) + s2c(i,j,4)*q0(i1,j1+1,k) enddo enddo @@ -1083,7 +1033,8 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) if (is_master()) call pmaxmin( 'MOIST_Q_', q0(:,:,:), im*jm, km, 1.0_FVPRC) enddo if (filetype == 0) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) else close(IUNIT) end if @@ -1124,21 +1075,20 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) if (filetype /= 0) then open(IUNIT,file=triM(tracer_bundles(ifile)%file_name) ,access='sequential',form='unformatted',status='old') else - NCIO = MAPL_NCIOOpen(trim(tracer_bundles(ifile)%file_name)) - call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars) + call formatter%open(trim(tracer_bundles(ifile)%file_name),pFIO_READ,rc=status) end if do ivar=1,size(tracer_bundles(ifile)%vars) nlev=tracer_bundles(ifile)%vars(ivar)%nLev - if (filetype == 0) call MAPL_NCIOGetVarName(NCIO,ivar,vname) + if (filetype == 0) vname = trim(tracer_bundles(ifile)%vars(ivar)%name) do k=1,nlev if (filetype /= 0) then read (IUNIT, IOSTAT=status)gslice_r4 else if (tracer_bundles(ifile)%vars(ivar)%nLev/=1) then - call MAPL_VarRead(NCIO,vname,gslice_r4,lev=k) + call MAPL_VarRead(formatter,vname,gslice_r4,lev=k) else - call MAPL_VarRead(NCIO,vname,gslice_r4) + call MAPL_VarRead(formatter,vname,gslice_r4) end if end if q0(1 :im/2,:,k) = gslice_r4(im/2 + 1 :im , :) ! Regrid from -180:180 to 0:360 @@ -1171,7 +1121,7 @@ subroutine get_geos_latlon_ic( Atm, extra_rst) enddo if (filetype == 0) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() deallocate(gslice_r4) end if @@ -1820,7 +1770,7 @@ end subroutine init_cubsph_grid subroutine interp_c2c_vect(npx_in, npy_in, npx_out, npy_out, npz, ntiles, domain_i, & is,ie, js,je, isd_i,ied_i, jsd_i,jed_i, is_i,ie_i, js_i,je_i, & - u_in, v_in, u_out, v_out, index_c2c, weight_c2c, corner_in, corner_out, gridstruct) + u_in, v_in, u_out, v_out, regridder, corner_in, corner_out, gridstruct) use GRID_UTILS_mod, only: latlon2xyz use GRID_UTILS_mod, only: get_dx, get_dxa, get_dy, get_dya, & get_center_vect, get_west_vect, & @@ -1831,8 +1781,7 @@ subroutine interp_c2c_vect(npx_in, npy_in, npx_out, npy_out, npz, ntiles, domain integer, intent(IN) :: npx_in, npy_in, npx_out, npy_out, npz, ntiles type(domain2d), intent(INOUT) :: domain_i integer, intent(IN) :: is,ie, js,je, isd_i,ied_i, jsd_i,jed_i, is_i,ie_i, js_i,je_i - integer, intent(IN) :: index_c2c(3, is:ie, js:je ) - real(REAL64), intent(IN) :: weight_c2c(4, is:ie, js:je ) + class(AbstractRegridder), pointer :: regridder real(REAL64), intent(IN) :: corner_in(2,is_i:ie_i+1,js_i:je_i+1) real(REAL64), intent(IN) :: corner_out(2,is:ie+1,js:je+1) real(FVPRC), intent(IN) :: u_in(isd_i:ied_i,jsd_i:jed_i+1,npz) @@ -1843,8 +1792,9 @@ subroutine interp_c2c_vect(npx_in, npy_in, npx_out, npy_out, npz, ntiles, domain real(FVPRC) :: tmp(isd_i:ied_i,jsd_i:jed_i) - integer :: i,j,l,k,n + integer :: i,j,k,n real(REAL64), dimension(:,:,:), allocatable :: vxyz_in, vxyz_out + real(REAL32), dimension(:,:), allocatable :: tmp_in, tmp_out real(REAL64), dimension(:,:,:), allocatable :: ec1, ec2, ew1, ew2, es1, es2 real(REAL64), dimension(:,:) , allocatable :: dx, dy, dxa, dya, rdxa, rdya, cosa_s, sina_s real(FVPRC) :: u1, v1, vx, vy, vz @@ -1882,6 +1832,8 @@ subroutine interp_c2c_vect(npx_in, npy_in, npx_out, npy_out, npz, ntiles, domain allocate(ec1(3,isd_i:ied_i,jsd_i:jed_i), ec2(3,isd_i:ied_i,jsd_i:jed_i)) allocate(cosa_s(isd_i:ied_i,jsd_i:jed_i), sina_s(isd_i:ied_i,jsd_i:jed_i)) allocate(vxyz_out(3,is:ie,js:je)) + allocate(tmp_in(is_i:ie_i,js_i:je_i)) + allocate(tmp_out(is:ie,js:je)) !-------------------------------------------------------! ! geometrical properties of input grid ! @@ -1924,19 +1876,9 @@ subroutine interp_c2c_vect(npx_in, npy_in, npx_out, npy_out, npz, ntiles, domain vxyz_in(n,:,:) = tmp enddo do n=1,3 - do j=js,je - do i=is,ie -!----------------------------------------------------------! -! do interpolation of flow vector on A-Grid ! -!----------------------------------------------------------! - ic=index_c2c(1,i,j) - jc=index_c2c(2,i,j) - vxyz_out(n,i,j)=weight_c2c(1,i,j)*vxyz_in(n,ic ,jc ) & - +weight_c2c(2,i,j)*vxyz_in(n,ic ,jc+1) & - +weight_c2c(3,i,j)*vxyz_in(n,ic+1,jc+1) & - +weight_c2c(4,i,j)*vxyz_in(n,ic+1,jc ) - enddo - enddo + tmp_in=vxyz_in(n,is_i:ie_i,js_i:je_i) + call regridder%regrid(tmp_in,tmp_out,rc=status) + vxyz_out(n,:,:)=tmp_out enddo do j=js,je do i=is,ie diff --git a/geos_utils/fv_regridding_utils.F90 b/geos_utils/fv_regridding_utils.F90 index 68086d494..5cb009141 100644 --- a/geos_utils/fv_regridding_utils.F90 +++ b/geos_utils/fv_regridding_utils.F90 @@ -8,10 +8,7 @@ module fv_regridding_utils use fv_mp_mod, only: is_master, ng use fv_mapz_mod, only: mappm use mpp_mod, only: mpp_error, FATAL, NOTE, mpp_broadcast,mpp_npes - !use MAPL_MOD, only: MAPL_PI_R8, MAPL_OMEGA, MAPL_GRAV, & - !MAPL_KAPPA, MAPL_RGAS, MAPL_RVAP, & - !MAPL_CP - use MAPL_MOD + use MAPL implicit none @@ -20,7 +17,6 @@ module fv_regridding_utils public remap_scalar public fv_rst public copy_fv_rst - public simple_cs_grid_creator real(FVPRC), parameter :: PI = MAPL_PI_R8 real(FVPRC), parameter :: OMEGA = MAPL_OMEGA @@ -34,8 +30,8 @@ module fv_regridding_utils type fv_var character(len=128) :: name integer :: nlev - real(REAL8), pointer :: ptr2d(:,:) => null() - real(REAL8), pointer :: ptr3d(:,:,:) => null() + real(FVPRC), pointer :: ptr2d(:,:) => null() + real(FVPRC), pointer :: ptr3d(:,:,:) => null() end type fv_var type fv_rst @@ -254,43 +250,5 @@ subroutine remap_edge(q1,q2,is,ie,km,kn,ak,bk) end subroutine remap_edge -function simple_cs_grid_creator(IM_World,JM_World,NX,NY,rc) result(esmfgrid) - - integer, intent(IN) :: IM_WORLD, JM_WORLD - integer, intent(IN) :: NX, NY - integer, optional, intent(OUT) :: rc - type (ESMF_Grid) :: esmfgrid - - integer, allocatable :: IMS(:), JMS(:) - integer :: n - integer :: STATUS - character(len=ESMF_MAXSTR), parameter :: Iam="simple_cs_grid_creator" - - allocate( IMS(0:NX-1) ) - allocate( JMS(0:NY-1) ) - - call MAPL_DecomposeDim ( IM_WORLD , IMS , NX , symmetric=.true. ) - call MAPL_DecomposeDim ( JM_WORLD/6, JMS(0:NY/6 -1) , NY/6, symmetric=.true. ) - do n=2,6 - JMS((n-1)*NY/6 : n*NY/6 -1) = JMS(0:NY/6 -1) - enddo - - esmfgrid = ESMF_GridCreate( & - name="dummy", & - countsPerDEDim1=IMS, & - countsPerDEDim2=JMS, & - indexFlag = ESMF_INDEX_USER, & - gridMemLBound = (/1,1/), & - gridEdgeLWidth = (/0,0/), & - gridEdgeUWidth = (/0,0/), & - coordDep1 = (/1,2/), & - coordDep2 = (/1,2/), & - rc=status) - VERIFY_(STATUS) - - RETURN_(ESMF_SUCCESS) - -end function simple_cs_grid_creator - end module fv_regridding_utils diff --git a/geos_utils/rs_scaleMod.F90 b/geos_utils/rs_scaleMod.F90 index 8f85a2b9b..22d5c6987 100644 --- a/geos_utils/rs_scaleMod.F90 +++ b/geos_utils/rs_scaleMod.F90 @@ -2,10 +2,9 @@ module rs_scaleMod use fv_arrays_mod - use MAPL_ConstantsMod, only: MAPL_PSDRY - use MAPL_Mod + use MAPL + use gFTL_StringIntegerMap use ESMF - use pFIO_StringIntegerMapMod use, intrinsic :: iso_fortran_env, only: REAL64, REAL32 ! bma added implicit none diff --git a/model/mapz-driver/CMakeLists.txt b/model/mapz-driver/CMakeLists.txt index f0e6de384..3aced14f5 100644 --- a/model/mapz-driver/CMakeLists.txt +++ b/model/mapz-driver/CMakeLists.txt @@ -33,7 +33,7 @@ endif () ecbuild_add_executable( TARGET mapz-driver SOURCES ${srcs} - LIBS ${FMS} ${MAPL_BASE}) + LIBS ${FMS} MAPL) target_compile_definitions(mapz-driver PRIVATE MAPL_MODE SPMD TIMING) target_compile_options(mapz-driver PRIVATE ${TRACEBACK}) set_target_properties(${this} PROPERTIES Fortran_MODULE_DIRECTORY ${esma_include}/${this}) diff --git a/tools/external_ic.F90 b/tools/external_ic.F90 index 4fedbe293..6c3966fac 100644 --- a/tools/external_ic.F90 +++ b/tools/external_ic.F90 @@ -1,7 +1,7 @@ module external_ic_mod #ifdef MAPL_MODE -#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 +#define DEALLOCGLOB_(A) if(associated(A)) then;A=0;if(MAPL_ShmInitialized) then; call MAPL_DeAllocNodeArray(A,rc=status);else; deallocate(A);endif;NULLIFY(A);endif #endif #ifndef DYCORE_SOLO @@ -18,8 +18,7 @@ module external_ic_mod use constants_mod, only: pi=>pi_8, omega, grav, kappa, rdgas, rvgas, cp_air #ifdef MAPL_MODE - use MAPL_IOMod - use MAPL_ShmemMod + use MAPL #endif use, intrinsic :: iso_fortran_env, only: REAL64, REAL32 @@ -482,12 +481,18 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) integer :: filetype logical :: isNC4 - type(MAPL_NCIO) :: ncio + type(Netcdf4_Fileformatter) :: formatter + type(FileMetadata), allocatable :: cfg(:) integer :: nDims, nVars, ivar, dimSizes(3) character(len=128) :: vname real(FVPRC), allocatable :: gslice_r4(:,:) real*8, allocatable :: gslice_r8(:,:) integer :: tileoff,lvar_cnt + character(len=128), allocatable :: vnames(:) + character(len=:), pointer :: var_name + type(StringVariableMap), pointer :: vars + type(StringVariableMapIterator) :: iter + !bma added character(len=128) :: moist_order(9) = (/"Q ","QLLS","QLCN","CLLS","CLCN","QILS","QICN","NCPL","NCPI"/) @@ -523,8 +528,13 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) if (isNC4) then - NCIO = MAPL_NCIOOpen(fname) - call MAPL_NCIOGetDimSizes(NCIO,lon=im,lat=jm,lev=km) + allocate(cfg(1)) + call formatter%open(fname,pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + im =cfg(1)%get_dimension('lon',rc=status) + jm =cfg(1)%get_dimension('lat',rc=status) + km =cfg(1)%get_dimension('lev',rc=status) + allocate(gslice_r8(im,jm)) else @@ -594,13 +604,13 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) allocate ( bk0(km+1) ) allocate ( akbk_r8(km+1) ) if (isNC4) then - call MAPL_VarRead(NCIO,"AK",akbk_r8) + call MAPL_VarRead(formatter,"AK",akbk_r8) else read (IUNIT, IOSTAT=status) akbk_r8 end if ak0 = akbk_r8 if (isNC4) then - call MAPL_VarRead(NCIO,"BK",akbk_r8) + call MAPL_VarRead(formatter,"BK",akbk_r8) else read (IUNIT, IOSTAT=status) akbk_r8 end if @@ -615,7 +625,7 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) if (isNC4) then tileoff = (tile-1)*(jm/ntiles) do k=1,km - call MAPL_VarRead(NCIO,"U",gslice_r8,lev=k) + call MAPL_VarRead(formatter,"U",gslice_r8,lev=k) u0(is_i:ie_i,js_i:je_i,k) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i) enddo else @@ -631,7 +641,7 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) if (isNC4) then tileoff = (tile-1)*(jm/ntiles) do k=1,km - call MAPL_VarRead(NCIO,"V",gslice_r8,lev=k) + call MAPL_VarRead(formatter,"V",gslice_r8,lev=k) v0(is_i:ie_i,js_i:je_i,k) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i) enddo else @@ -677,7 +687,7 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) if (isNC4) then tileoff = (tile-1)*(jm/ntiles) do k=1,km - call MAPL_VarRead(NCIO,"PT",gslice_r8,lev=k) + call MAPL_VarRead(formatter,"PT",gslice_r8,lev=k) t0(is_i:ie_i,js_i:je_i,k) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i) enddo else @@ -690,7 +700,7 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) ps0(:,:) = 0.0 if (isNC4) then tileoff = (tile-1)*(jm/ntiles) - call MAPL_VarRead(NCIO,"PE",gslice_r8,lev=km+1) + call MAPL_VarRead(formatter,"PE",gslice_r8,lev=km+1) ps0(is_i:ie_i,js_i:je_i) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i) else slice_2d = npts*npts*ntiles @@ -707,7 +717,7 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) if (isNC4) then tileoff = (tile-1)*(jm/ntiles) do k=1,km - call MAPL_VarRead(NCIO,"PKZ",gslice_r8,lev=k) + call MAPL_VarRead(formatter,"PKZ",gslice_r8,lev=k) pkz0(is_i:ie_i,js_i:je_i) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i) t0(is_i:ie_i,js_i:je_i,k) = t0(is_i:ie_i,js_i:je_i,k)*pkz0(is_i:ie_i,js_i:je_i) enddo @@ -723,7 +733,8 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) deallocate ( pkz0 ) if (isNC4) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) deallocate(gslice_r8) end if @@ -798,8 +809,10 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) else lvar_cnt = 0 allocate(gslice_r4(im,jm)) - NCIO = MAPL_NCIOOpen("moist_internal_restart_in") - call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars) + allocate(cfg(1)) + call formatter%open("moist_internal_restart_in",pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + call MAPL_IOCountNonDimVars(cfg(1),nvars,rc=status) if (nVars /= iq_moist1-iq_moist0+1) call mpp_error(FATAL,'Wrong number of variables in moist file') tileoff = (tile-1)*(jm/ntiles) end if @@ -812,7 +825,7 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) call mpp_update_domains(q0(:,:,k), domain_i) else vname = trim(moist_order(lvar_cnt)) - call MAPL_VarRead(NCIO,vname,gslice_r4,lev=k) + call MAPL_VarRead(formatter,vname,gslice_r4,lev=k) q0(is_i:ie_i,js_i:je_i,k)=gslice_r4(is_i:ie+i,tileoff+js_i:tileoff+je_i) end if call mpp_update_domains(q0(:,:,k), domain_i) @@ -833,7 +846,8 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) enddo if (filetype == 0) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) deallocate(gslice_r4) end if @@ -852,12 +866,29 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) else lvar_cnt = 0 allocate(gslice_r4(im,jm)) - NCIO = MAPL_NCIOOpen("gocart_internal_restart_in") - call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars) + allocate(cfg(1)) + call formatter%open("gocart_internal_restart_in",pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + call MAPL_IOCountNonDimVars(cfg(1),nvars,rc=status) if (nVars /= iq_gocart1-iq_gocart0+1) call mpp_error(FATAL,'Wrong number of variables in gocart file') tileoff = (tile-1)*(jm/ntiles) + + allocate(vnames(nVars)) + vars => cfg(1)%get_variables() + iter = vars%begin() + + lvar_cnt=0 + do while(iter /= vars%end()) + var_name => iter%key() + if (.not.cfg(1)%is_coordinate_variable(var_name)) then + lvar_cnt=lvar_cnt+1 + vnames(lvar_cnt)=var_name + end if + call iter%next() + enddo end if + lvar_cnt = 0 do ivar=iq_gocart0,iq_gocart1 if (filetype ==0) lvar_cnt=lvar_cnt+1 do k=1,km @@ -865,8 +896,7 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) call parallel_read_file_r4('gocart_internal_restart_in', npts, is_i,ie_i, js_i,je_i, 1, offset, q0(is_i:ie_i,js_i:je_i,k)) call mpp_update_domains(q0(:,:,k), domain_i) else - call MAPL_NCIOGetVarName(NCIO,lvar_cnt,vname) - call MAPL_VarRead(NCIO,vname,gslice_r4,lev=k) + call MAPL_VarRead(formatter,vnames(lvar_cnt),gslice_r4,lev=k) q0(is_i:ie_i,js_i:je_i,k)=gslice_r4(is_i:ie+i,tileoff+js_i:tileoff+je_i) end if call mpp_update_domains(q0(:,:,k), domain_i) @@ -886,7 +916,9 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) enddo if (filetype == 0) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) + deallocate(vnames) deallocate(gslice_r4) end if @@ -904,12 +936,29 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) else lvar_cnt = 0 allocate(gslice_r4(im,jm)) - NCIO = MAPL_NCIOOpen("pchem_internal_restart_in") - call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars) + allocate(cfg(1)) + call formatter%open("pchem_internal_restart_in",pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + call MAPL_IOCountNonDimVars(cfg(1),nvars,rc=status) if (nVars /= iq_pchem1-iq_pchem0+1) call mpp_error(FATAL,'Wrong number of variables in pchem file') tileoff = (tile-1)*(jm/ntiles) + + allocate(vnames(nVars)) + vars => cfg(1)%get_variables() + iter = vars%begin() + + lvar_cnt=0 + do while(iter /= vars%end()) + var_name => iter%key() + if (.not.cfg(1)%is_coordinate_variable(var_name)) then + lvar_cnt=lvar_cnt+1 + vnames(lvar_cnt)=var_name + end if + call iter%next() + enddo end if + lvar_cnt=0 do ivar=iq_pchem0,iq_pchem1 if (filetype == 0) lvar_cnt=lvar_cnt+1 do k=1,km @@ -917,8 +966,7 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) call parallel_read_file_r4('pchem_internal_restart_in', npts, is_i,ie_i, js_i,je_i, 1, offset, q0(is_i:ie_i,js_i:je_i,k)) call mpp_update_domains(q0(:,:,k), domain_i) else - call MAPL_NCIOGetVarName(NCIO,lvar_cnt,vname) - call MAPL_VarRead(NCIO,vname,gslice_r4,lev=k) + call MAPL_VarRead(formatter,vnames(lvar_cnt),gslice_r4,lev=k) q0(is_i:ie_i,js_i:je_i,k)=gslice_r4(is_i:ie+i,tileoff+js_i:tileoff+je_i) end if call mpp_update_domains(q0(:,:,k), domain_i) @@ -938,7 +986,9 @@ subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers ) enddo if (filetype == 0) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) + deallocate(vnames) deallocate(gslice_r4) end if @@ -1046,13 +1096,20 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) integer :: filetype logical :: isNC4 - type(MAPL_NCIO) :: ncio integer :: nDims, nVars, ivar, dimSizes(3) character(len=128) :: vname integer :: iq_moist0 , iq_moist1 integer :: iq_gocart0, iq_gocart1 integer :: iq_pchem0 , iq_pchem1 integer :: lvar_cnt + type(Netcdf4_Fileformatter) :: formatter + type(FileMetadata), allocatable :: cfg(:) + character(len=128), allocatable :: vnames(:) + character(len=:), pointer :: var_name + type(StringVariableMap), pointer :: vars + type(StringVariableMapIterator) :: iter + + !bma added character(len=128) :: moist_order(9) = (/"Q ","QLLS","QLCN","CLLS","CLCN","QILS","QICN","NCPL","NCPI"/) @@ -1086,8 +1143,12 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) if (isNC4) then - NCIO = MAPL_NCIOOpen(fname) - call MAPL_NCIOGetDimSizes(NCIO,lon=im,lat=jm,lev=km) + allocate(cfg(1)) + call formatter%open(fname,pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + im =cfg(1)%get_dimension('lon',rc=status) + jm =cfg(1)%get_dimension('lat',rc=status) + km =cfg(1)%get_dimension('lev',rc=status) else @@ -1121,13 +1182,13 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) allocate ( bk0(km+1) ) allocate ( akbk_r8(km+1) ) if (isNC4) then - call MAPL_VarRead(NCIO,"AK",akbk_r8) + call MAPL_VarRead(formatter,"AK",akbk_r8) else read (IUNIT, IOSTAT=status) akbk_r8 end if ak0 = akbk_r8 if (isNC4) then - call MAPL_VarRead(NCIO,"BK",akbk_r8) + call MAPL_VarRead(formatter,"BK",akbk_r8) else read (IUNIT, IOSTAT=status) akbk_r8 end if @@ -1141,7 +1202,7 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) allocate ( u0(im,jm,km) ) do k=1,km if (isNC4) then - call MAPL_VarRead(NCIO,"U",r8latlon,lev=k) + call MAPL_VarRead(formatter,"U",r8latlon,lev=k) else read (IUNIT, IOSTAT=status) r8latlon end if @@ -1154,7 +1215,7 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) allocate ( v0(im,jm,km) ) do k=1,km if (isNC4) then - call MAPL_VarRead(NCIO,"V",r8latlon,lev=k) + call MAPL_VarRead(formatter,"V",r8latlon,lev=k) else read (IUNIT, IOSTAT=status) r8latlon end if @@ -1202,7 +1263,7 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) allocate ( t0(im,jm,km) ) do k=1,km if (isNC4) then - call MAPL_VarRead(NCIO,"PT",r8latlon,lev=k) + call MAPL_VarRead(formatter,"PT",r8latlon,lev=k) else read (IUNIT, IOSTAT=status) r8latlon end if @@ -1214,7 +1275,7 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) ! Read PE do k=1,km+1 if (isNC4) then - call MAPL_VarRead(NCIO,"PE",r8latlon,lev=k) + call MAPL_VarRead(formatter,"PE",r8latlon,lev=k) else read (IUNIT, IOSTAT=status) r8latlon end if @@ -1226,7 +1287,7 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) allocate ( pkz0(im,jm) ) do k=1,km if (isNC4) then - call MAPL_VarRead(NCIO,"PKZ",r8latlon,lev=k) + call MAPL_VarRead(formatter,"PKZ",r8latlon,lev=k) else read (IUNIT, IOSTAT=status) r8latlon end if @@ -1240,7 +1301,8 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) call print_memuse_stats('get_geos_latlon_ic: converted T') deallocate ( pkz0 ) if (isNC4) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) else close (IUNIT) end if @@ -1329,8 +1391,9 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) open(IUNIT,file="moist_internal_restart_in" ,access='sequential',form='unformatted',status='old') else lvar_cnt = 0 - NCIO = MAPL_NCIOOpen("moist_internal_restart_in") - call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars) + call formatter%open("moist_internal_restart_in",pFIO_READ,rc=status) + cfg = formatter%read(rc=status) + call MAPL_IOCountNonDimVars(cfg(1),nvars,rc=status) if (nVars /= iq_moist1-iq_moist0+1) call mpp_error(FATAL,'Wrong number of variables in moist file') end if @@ -1341,7 +1404,7 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) read (IUNIT, IOSTAT=status) r4latlon else vname = trim(moist_order(lvar_cnt)) - call MAPL_VarRead(NCIO,vname,r4latlon,lev=k) + call MAPL_VarRead(formatter,vname,r4latlon,lev=k) end if q0(1 :im/2,:,k) = r4latlon(im/2 + 1 :im , :) ! Regrid from -180:180 to 0:360 q0(im/2 + 1:im ,:,k) = r4latlon(1 :im/2, :) ! Regrid from -180:180 to 0:360 @@ -1359,7 +1422,8 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) if (is_master()) call pmaxmin( 'MOIST_Q_', q0(:,:,:), im*jm, km, 1.0_FVPRC) enddo if (filetype == 0) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) else close(IUNIT) end if @@ -1379,19 +1443,34 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) open(IUNIT,file="gocart_internal_restart_in" ,access='sequential',form='unformatted',status='old') else lvar_cnt = 0 - NCIO = MAPL_NCIOOpen("gocart_internal_restart_in") - call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars) + allocate(cfg(1)) + call formatter%open("gocart_internal_restart_in",pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + call MAPL_IOCountNonDimVars(cfg(1),nvars,rc=status) if (nVars /= iq_gocart1-iq_gocart0+1) call mpp_error(FATAL,'Wrong number of variables in gocart file') + + allocate(vnames(nVars)) + vars => cfg(1)%get_variables() + iter = vars%begin() + + do while(iter /= vars%end()) + var_name => iter%key() + if (.not.cfg(1)%is_coordinate_variable(var_name)) then + lvar_cnt=lvar_cnt+1 + vnames(lvar_cnt)=var_name + end if + enddo + end if + lvar_cnt=0 do ivar=iq_gocart0,iq_gocart1 if (filetype ==0) lvar_cnt=lvar_cnt+1 do k=1,km if (filetype /= 0) then read (IUNIT, IOSTAT=status) r4latlon else - call MAPL_NCIOGetVarName(NCIO,lvar_cnt,vname) - call MAPL_VarRead(NCIO,vname,r4latlon,lev=k) + call MAPL_VarRead(formatter,vnames(lvar_cnt),r4latlon,lev=k) end if q0(1 :im/2,:,k) = r4latlon(im/2 + 1 :im , :) ! Regrid from -180:180 to 0:360 q0(im/2 + 1:im ,:,k) = r4latlon(1 :im/2, :) ! Regrid from -180:180 to 0:360 @@ -1410,7 +1489,9 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) enddo if (filetype == 0) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) + deallocate(vnames) else close(IUNIT) end if @@ -1429,9 +1510,24 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) open(IUNIT,file="pchem_internal_restart_in" ,access='sequential',form='unformatted',status='old') else lvar_cnt = 0 - NCIO = MAPL_NCIOOpen("pchem_internal_restart_in") - call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars) + allocate(cfg(1)) + call formatter%open("pchem_internal_restart_in",pFIO_READ,rc=status) + cfg(1) = formatter%read(rc=status) + call MAPL_IOCountNonDimVars(cfg(1),nvars,rc=status) if (nVars /= iq_pchem1-iq_pchem0+1) call mpp_error(FATAL,'Wrong number of variables in pchem file') + + allocate(vnames(nVars)) + vars => cfg(1)%get_variables() + iter = vars%begin() + + do while(iter /= vars%end()) + var_name => iter%key() + if (.not.cfg(1)%is_coordinate_variable(var_name)) then + lvar_cnt=lvar_cnt+1 + vnames(lvar_cnt)=var_name + end if + enddo + end if do ivar=iq_pchem0,iq_pchem1 @@ -1440,8 +1536,7 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) if (filetype /= 0) then read (IUNIT, IOSTAT=status) r4latlon else - call MAPL_NCIOGetVarName(NCIO,lvar_cnt,vname) - call MAPL_VarRead(NCIO,vname,r4latlon,lev=k) + call MAPL_VarRead(formatter,vnames(lvar_cnt),r4latlon,lev=k) end if q0(1 :im/2,:,k) = r4latlon(im/2 + 1 :im , :) ! Regrid from -180:180 to 0:360 q0(im/2 + 1:im ,:,k) = r4latlon(1 :im/2, :) ! Regrid from -180:180 to 0:360 @@ -1460,7 +1555,9 @@ subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers) enddo if (filetype == 0) then - call MAPL_NCIOClose(NCIO,destroy=.true.) + call formatter%close() + deallocate(cfg) + deallocate(vnames) else close(IUNIT) end if @@ -2827,6 +2924,7 @@ subroutine interp_c2c_vect(npx_in, npy_in, npx_out, npy_out, npz, ntiles, domain !------------------------------------------------------------------! ! calculate xyz cell corners and cell centers ! !------------------------------------------------------------------! + tmp = 0.0 allocate(xyz_corner_in (3, isd_i:ied_i+1, jsd_i:jed_i+1), & xyz_corner_out(3, is :ie +1, js :je +1)) do j=js_i,je_i+1 diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index df42311fc..f53e93404 100644 --- a/tools/fv_treat_da_inc.F90 +++ b/tools/fv_treat_da_inc.F90 @@ -119,7 +119,7 @@ module fv_treat_da_inc_mod use field_manager_mod, only: MODEL_ATMOS #ifdef MAPL_MODE - use MAPL_MOD + use MAPL #else use constants_mod, only: pi=>pi_8, omega, grav, kappa, & rdgas, rvgas, cp_air diff --git a/tools/init_hydro.F90 b/tools/init_hydro.F90 index bbdc42a5c..14288336f 100644 --- a/tools/init_hydro.F90 +++ b/tools/init_hydro.F90 @@ -173,25 +173,23 @@ subroutine p_var(km, ifirst, ilast, jfirst, jlast, ptop, ptop_min, & enddo - if ( .not.hydrostatic ) then - - rdg = -rdgas / grav - if ( present(make_nh) ) then - if ( make_nh ) then - delz = 1.e25 -!$OMP parallel do default(none) shared(ifirst,ilast,jfirst,jlast,km,delz,rdg,pt,peln) - do k=1,km - do j=jfirst,jlast - do i=ifirst,ilast - delz(i,j,k) = rdg*pt(i,j,k)*(peln(i,k+1,j)-peln(i,k,j)) - enddo - enddo - enddo - if(is_master()) write(*,*) 'delz computed from hydrostatic state' - endif - endif - -! This is a bug below. Clearly this needs to be protected as we can not use delz if it was not computed above + !if ( .not.hydrostatic ) then + + !rdg = -rdgas / grav + !if ( present(make_nh) ) then + !if ( make_nh ) then + !delz = 1.e25 +!!$OMP parallel do default(none) shared(ifirst,ilast,jfirst,jlast,km,delz,rdg,pt,peln) + !do k=1,km + !do j=jfirst,jlast + !do i=ifirst,ilast + !delz(i,j,k) = rdg*pt(i,j,k)*(peln(i,k+1,j)-peln(i,k,j)) + !enddo + !enddo + !enddo + !if(is_master()) write(*,*) 'delz computed from hydrostatic state' + !endif + !endif !if ( moist_phys ) then !!------------------------------------------------------------------ @@ -225,7 +223,7 @@ subroutine p_var(km, ifirst, ilast, jfirst, jlast, ptop, ptop_min, & !enddo !endif - endif + !endif end subroutine p_var