diff --git a/mediator/med_io_mod.F90 b/mediator/med_io_mod.F90 index c86f87c7..2365dce4 100644 --- a/mediator/med_io_mod.F90 +++ b/mediator/med_io_mod.F90 @@ -729,7 +729,7 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, & type(ESMF_CoordSys_Flag) :: coordsys integer :: rcode integer :: nf,ns,ng - integer :: k,n + integer :: k,n,n2 integer :: ndims, nelements integer ,target :: dimid2(2) integer ,target :: dimid3(3) @@ -754,15 +754,22 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, & integer, pointer :: Dof(:) real(r8), pointer :: fldptr1(:) real(r8), pointer :: fldptr2(:,:) + real(r8), pointer :: fldptr3(:,:,:) real(r8), allocatable :: ownedElemCoords(:), ownedElemCoords_x(:), ownedElemCoords_y(:) character(CS) :: cnumber + character(CS) :: cnumber2 character(CL) :: tmpstr type(ESMF_Field) :: lfield integer :: rank - integer :: ungriddedUBound(1) ! currently the size must equal 1 for rank 2 fields - integer :: gridToFieldMap(1) ! currently the size must equal 1 for rank 2 fields logical :: tiles character(CL), allocatable :: fieldNameList(:) + + ! For a single ungridded dimension, there will be 1 element in ungriddedUBound and 1 + ! element in gridToFieldMap; for two ungridded dimensions, there will be 2 elements in + ! ungriddedUBound but still 1 element in gridToFieldMap. + integer :: ungriddedUBound(2) + integer :: gridToFieldMap(1) + character(*),parameter :: subName = '(med_io_write_FB) ' !------------------------------------------------------------------------------- @@ -935,12 +942,47 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, & ! TODO (mvertens, 2019-03-13): this is a temporary mod to NOT write hgt if (trim(itemc) /= "hgt") then - if (rank == 2) then + if (rank == 3) then call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return write(cnumber,'(i0)') ungriddedUbound(1) + write(cnumber2,'(i0)') ungriddedUbound(2) call ESMF_LogWrite(trim(subname)//':'//'field '//trim(itemc)// & - ' has an griddedUBound of '//trim(cnumber), ESMF_LOGMSG_INFO) + ' has an ungriddedUBound of '//trim(cnumber)//' x '//trim(cnumber2), ESMF_LOGMSG_INFO) + + ! Create a new output variable for each element of the 2 ungridded dimensions + do n2 = 1,ungriddedUBound(2) + do n = 1,ungriddedUBound(1) + write(cnumber,'(i0)') n + write(cnumber2,'(i0)') n2 + name1 = trim(lpre)//'_'//trim(itemc)//trim(cnumber)//'_'//trim(cnumber2) + call ESMF_LogWrite(trim(subname)//': defining '//trim(name1), ESMF_LOGMSG_INFO) + if (luse_float) then + rcode = pio_def_var(io_file, trim(name1), PIO_REAL, dimid, varid) + rcode = pio_put_att(io_file, varid,"_FillValue",real(lfillvalue,r4)) + else + rcode = pio_def_var(io_file, trim(name1), PIO_DOUBLE, dimid, varid) + rcode = pio_put_att(io_file,varid,"_FillValue",lfillvalue) + end if + if (NUOPC_FieldDictionaryHasEntry(trim(itemc))) then + call NUOPC_FieldDictionaryGetEntry(itemc, canonicalUnits=cunit, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + rcode = pio_put_att(io_file, varid, "units", trim(cunit)) + end if + rcode = pio_put_att(io_file, varid, "standard_name", trim(name1)) + if (present(tavg)) then + if (tavg) then + rcode = pio_put_att(io_file, varid, "cell_methods", "time: mean") + endif + endif + end do + end do + else if (rank == 2) then + call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + write(cnumber,'(i0)') ungriddedUbound(1) + call ESMF_LogWrite(trim(subname)//':'//'field '//trim(itemc)// & + ' has an ungriddedUBound of '//trim(cnumber), ESMF_LOGMSG_INFO) ! Create a new output variable for each element of the undistributed dimension do n = 1,ungriddedUBound(1) @@ -958,7 +1000,7 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, & if (NUOPC_FieldDictionaryHasEntry(trim(itemc))) then call NUOPC_FieldDictionaryGetEntry(itemc, canonicalUnits=cunit, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - rcode = pio_put_att(io_file, varid, "units" , trim(cunit)) + rcode = pio_put_att(io_file, varid, "units", trim(cunit)) end if rcode = pio_put_att(io_file, varid, "standard_name", trim(name1)) if (present(tavg)) then @@ -968,15 +1010,15 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, & endif end if end do - else + else if (rank == 1) then name1 = trim(lpre)//'_'//trim(itemc) - call ESMF_LogWrite(trim(subname)//':'//trim(itemc)//':'//trim(name1),ESMF_LOGMSG_INFO) + call ESMF_LogWrite(trim(subname)//': defining '//trim(name1), ESMF_LOGMSG_INFO) if (luse_float) then rcode = pio_def_var(io_file, trim(name1), PIO_REAL, dimid, varid) - rcode = pio_put_att(io_file, varid, "_FillValue", real(lfillvalue, r4)) + rcode = pio_put_att(io_file, varid,"_FillValue",real(lfillvalue,r4)) else rcode = pio_def_var(io_file, trim(name1), PIO_DOUBLE, dimid, varid) - rcode = pio_put_att(io_file, varid, "_FillValue", lfillvalue) + rcode = pio_put_att(io_file,varid,"_FillValue",lfillvalue) end if if (NUOPC_FieldDictionaryHasEntry(trim(itemc))) then call NUOPC_FieldDictionaryGetEntry(itemc, canonicalUnits=cunit, rc=rc) @@ -988,7 +1030,10 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, & if (tavg) then rcode = pio_put_att(io_file, varid, "cell_methods", "time: mean") endif - end if + endif + else + call shr_log_error(subname//' ERROR: unhandled rank', line=__LINE__, file=u_FILE_u, rc=rc) + return end if end if end do @@ -1039,12 +1084,43 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, & end if call FB_getFldPtr(FB, itemc, & - fldptr1=fldptr1, fldptr2=fldptr2, rank=rank, rc=rc) + fldptr1=fldptr1, fldptr2=fldptr2, fldptr3=fldptr3, rank=rank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! TODO (mvertens, 2019-03-13): this is a temporary mod to NOT write hgt if (trim(itemc) /= "hgt") then - if (rank == 2) then + if (rank == 3) then + + ! Determine the size of the ungridded dimensions and the index where the distributed dimension is located + call ESMF_FieldBundleGet(FB, itemc, field=lfield, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, gridToFieldMap=gridToFieldMap, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (gridToFieldMap(1) /= 3) then + call shr_log_error( & + subname//' ERROR: For rank-3 fields, currently only gridToFieldMap(1)==3 is supported', & + line=__LINE__, file=u_FILE_u, rc=rc) + return + end if + + ! Output for each combination of ungriddedUbound indices + do n2 = 1,ungriddedUBound(2) + do n = 1,ungriddedUBound(1) + write(cnumber,'(i0)') n + write(cnumber2,'(i0)') n2 + name1 = trim(lpre)//'_'//trim(itemc)//trim(cnumber)//'_'//trim(cnumber2) + rcode = pio_inq_varid(io_file, trim(name1), varid) + call pio_setframe(io_file,varid,frame) + + if (luse_float) then + call pio_write_darray(io_file, varid, iodesc, real(fldptr3(n,n2,:),r4), rcode, fillval=real(lfillvalue,r4)) + else + call pio_write_darray(io_file, varid, iodesc, fldptr3(n,n2,:), rcode, fillval=lfillvalue) + end if + end do + end do + else if (rank == 2) then ! Determine the size of the ungridded dimension and the index where the undistributed dimension is located call ESMF_FieldBundleGet(FB, itemc, field=lfield, rc=rc) @@ -1489,7 +1565,7 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) type(ESMF_Field) :: lfield integer :: rcode integer :: nf - integer :: k,n,l + integer :: k,n,n2,l type(file_desc_t) :: pioid type(var_desc_t) :: varid type(io_desc_t) :: iodesc @@ -1500,11 +1576,18 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) integer :: rank, lsize real(r8), pointer :: fldptr1(:), fldptr1_tmp(:) real(r8), pointer :: fldptr2(:,:) + real(r8), pointer :: fldptr3(:,:,:) character(CL) :: tmpstr character(len=16) :: cnumber + character(len=16) :: cnumber2 integer(kind=Pio_Offset_Kind) :: lframe - integer :: ungriddedUBound(1) ! currently the size must equal 1 for rank 2 fieldds - integer :: gridToFieldMap(1) ! currently the size must equal 1 for rank 2 fieldds + + ! For a single ungridded dimension, there will be 1 element in ungriddedUBound and 1 + ! element in gridToFieldMap; for two ungridded dimensions, there will be 2 elements in + ! ungriddedUBound but still 1 element in gridToFieldMap. + integer :: ungriddedUBound(2) + integer :: gridToFieldMap(1) + character(*),parameter :: subName = '(med_io_read_FB) ' !------------------------------------------------------------------------------- rc = ESMF_Success @@ -1569,7 +1652,9 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_FieldGet(lfield, rank=rank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - if (rank == 2) then + if (rank == 3) then + name1 = trim(lpre)//'_'//trim(itemc)//'1_1' + else if (rank == 2) then name1 = trim(lpre)//'_'//trim(itemc)//'1' else if (rank == 1) then name1 = trim(lpre)//'_'//trim(itemc) @@ -1582,12 +1667,61 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! Get pointer to field bundle field - ! Field bundle might be 2d or 1d - but field on mediator history or restart file will always be 1d + ! Field bundle might be 3d, 2d or 1d - but field on mediator history or restart file will always be 1d call FB_getFldPtr(FB, itemc, & - fldptr1=fldptr1, fldptr2=fldptr2, rank=rank, rc=rc) + fldptr1=fldptr1, fldptr2=fldptr2, fldptr3=fldptr3, rank=rank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - if (rank == 2) then + if (rank == 3) then + + ! Determine the size of the ungridded dimensions and the + ! index where the distributed dimension is located + call ESMF_FieldBundleGet(FB, itemc, field=lfield, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, gridToFieldMap=gridToFieldMap, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (gridToFieldMap(1) /= 3) then + call shr_log_error( & + subname//' ERROR: For rank-3 fields, currently only gridToFieldMap(1)==3 is supported', & + line=__LINE__, file=u_FILE_u, rc=rc) + return + end if + + lsize = size(fldptr3, dim=3) + allocate(fldptr1_tmp(lsize)) + + do n2 = 1,ungriddedUBound(2) + do n = 1,ungriddedUBound(1) + ! Create a name for the 1d field on the mediator history or restart file based on the + ! ungridded dimension indices of the field bundle 3d field + write(cnumber,'(i0)') n + write(cnumber2,'(i0)') n2 + name1 = trim(lpre)//'_'//trim(itemc)//trim(cnumber)//'_'//trim(cnumber2) + + rcode = pio_inq_varid(pioid, trim(name1), varid) + if (rcode == pio_noerr) then + call ESMF_LogWrite(trim(subname)//' read field '//trim(name1), ESMF_LOGMSG_INFO) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call pio_setframe(pioid, varid, lframe) + call pio_read_darray(pioid, varid, iodesc, fldptr1_tmp, rcode) + rcode = pio_get_att(pioid, varid, "_FillValue", lfillvalue) + if (rcode /= pio_noerr) then + lfillvalue = fillvalue + endif + do l = 1,size(fldptr1_tmp) + if (fldptr1_tmp(l) == lfillvalue) fldptr1_tmp(l) = 0.0_r8 + enddo + else + fldptr1_tmp = 0.0_r8 + endif + fldptr3(n,n2,:) = fldptr1_tmp(:) + end do + end do + + deallocate(fldptr1_tmp) + + else if (rank == 2) then ! Determine the size of the ungridded dimension and the ! index where the undistributed dimension is located diff --git a/mediator/med_map_mod.F90 b/mediator/med_map_mod.F90 index a3d50c5e..2d5374c6 100644 --- a/mediator/med_map_mod.F90 +++ b/mediator/med_map_mod.F90 @@ -1300,6 +1300,9 @@ subroutine med_map_field_normalized(field_src, field_dst, routehandles, maptype, ! local variables integer :: n + real(r8), pointer :: data_src3d(:,:,:) + real(r8), pointer :: data_dst3d(:,:,:) + real(r8), pointer :: data_srctmp3d(:,:,:) real(r8), pointer :: data_src2d(:,:) real(r8), pointer :: data_dst2d(:,:) real(r8), pointer :: data_srctmp2d(:,:) @@ -1308,7 +1311,7 @@ subroutine med_map_field_normalized(field_src, field_dst, routehandles, maptype, real(r8), pointer :: data_srctmp1d(:) real(r8), pointer :: data_normsrc(:) real(r8), pointer :: data_normdst(:) - integer :: ungriddedUBound(1) ! currently the size must equal 1 for rank 2 fields + integer :: fieldrank integer :: lsize_src integer :: lsize_dst character(len=*), parameter :: subname=' (med_map_mod:med_map_field_normalized) ' @@ -1316,19 +1319,24 @@ subroutine med_map_field_normalized(field_src, field_dst, routehandles, maptype, rc = ESMF_SUCCESS - ! get a pointer (data_fracsrc) to the normalization array - ! get a pointer (data_src) to source field data in FBSrc - ! copy data_src to data_srctmp - - ! normalize data_src by data_fracsrc - call ESMF_FieldGet(field_normsrc, farrayPtr=data_normsrc, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return lsize_src = size(data_normsrc) - call ESMF_FieldGet(field_src, ungriddedUBound=ungriddedUBound, rc=rc) + call ESMF_FieldGet(field_src, rank=fieldrank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - if (ungriddedUbound(1) > 0) then + + ! Pre-normalization: multiply source data by the normalization field, + ! saving a copy so the source can be restored after regridding. + if (fieldrank == 3) then + call ESMF_FieldGet(field_src, farrayPtr=data_src3d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + allocate(data_srctmp3d(size(data_src3d,dim=1), size(data_src3d,dim=2), lsize_src)) + data_srctmp3d(:,:,:) = data_src3d(:,:,:) + do n = 1,lsize_src + data_src3d(:,:,n) = data_src3d(:,:,n) * data_normsrc(n) + end do + else if (fieldrank == 2) then call ESMF_FieldGet(field_src, farrayPtr=data_src2d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return allocate(data_srctmp2d(size(data_src2d,dim=1), lsize_src)) @@ -1336,7 +1344,7 @@ subroutine med_map_field_normalized(field_src, field_dst, routehandles, maptype, do n = 1,lsize_src data_src2d(:,n) = data_src2d(:,n) * data_normsrc(n) end do - else + else if (fieldrank == 1) then call ESMF_FieldGet(field_src, farrayPtr=data_src1d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return allocate(data_srctmp1d(lsize_src)) @@ -1344,32 +1352,52 @@ subroutine med_map_field_normalized(field_src, field_dst, routehandles, maptype, do n = 1,lsize_src data_src1d(n) = data_src1d(n) * data_normsrc(n) end do + else + call shr_log_error(subname//' ERROR: only ranks 1, 2, and 3 are supported', & + line=__LINE__, file=u_FILE_u, rc=rc) + return end if - ! regrid normalized packed source field + ! Regrid normalized source field to destination call med_map_field (field_src=field_src, field_dst=field_dst, routehandles=routehandles, maptype=maptype, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - ! restore original value to packed source field - if (ungriddedUbound(1) > 0) then + ! Restore original values to source field + if (fieldrank == 3) then + data_src3d(:,:,:) = data_srctmp3d(:,:,:) + deallocate(data_srctmp3d) + else if (fieldrank == 2) then data_src2d(:,:) = data_srctmp2d(:,:) deallocate(data_srctmp2d) - else + else if (fieldrank == 1) then data_src1d(:) = data_srctmp1d(:) deallocate(data_srctmp1d) + else + call shr_log_error(subname//' ERROR: only ranks 1, 2, and 3 are supported', & + line=__LINE__, file=u_FILE_u, rc=rc) + return end if - ! regrid normalization field from source to destination + ! Regrid normalization field from source to destination call med_map_field(field_src=field_normsrc, field_dst=field_normdst, routehandles=routehandles, maptype=maptype, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - ! get pointer to mapped fraction and normalize - ! destination mapped values by the reciprocal of the mapped fraction + ! Post-normalization: divide destination data by the mapped normalization field call ESMF_FieldGet(field_normdst, farrayPtr=data_normdst, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return lsize_dst = size(data_normdst) - if (ungriddedUbound(1) > 0) then + if (fieldrank == 3) then + call ESMF_FieldGet(field_dst, farrayPtr=data_dst3d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + do n = 1,lsize_dst + if (data_normdst(n) == 0.0_r8) then + data_dst3d(:,:,n) = 0.0_r8 + else + data_dst3d(:,:,n) = data_dst3d(:,:,n)/data_normdst(n) + end if + end do + else if (fieldrank == 2) then call ESMF_FieldGet(field_dst, farrayPtr=data_dst2d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return do n = 1,lsize_dst @@ -1379,7 +1407,7 @@ subroutine med_map_field_normalized(field_src, field_dst, routehandles, maptype, data_dst2d(:,n) = data_dst2d(:,n)/data_normdst(n) end if end do - else + else if (fieldrank == 1) then call ESMF_FieldGet(field_dst, farrayPtr=data_dst1d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return do n = 1,lsize_dst @@ -1389,6 +1417,10 @@ subroutine med_map_field_normalized(field_src, field_dst, routehandles, maptype, data_dst1d(n) = data_dst1d(n)/data_normdst(n) end if end do + else + call shr_log_error(subname//' ERROR: only ranks 1, 2, and 3 are supported', & + line=__LINE__, file=u_FILE_u, rc=rc) + return end if end subroutine med_map_field_normalized diff --git a/mediator/med_methods_mod.F90 b/mediator/med_methods_mod.F90 index b577b957..5ab81a49 100644 --- a/mediator/med_methods_mod.F90 +++ b/mediator/med_methods_mod.F90 @@ -22,12 +22,14 @@ module med_methods_mod interface med_methods_FieldPtr_compare ; module procedure & med_methods_FieldPtr_compare1, & - med_methods_FieldPtr_compare2 + med_methods_FieldPtr_compare2, & + med_methods_FieldPtr_compare3 end interface interface med_methods_check_for_nans module procedure med_methods_check_for_nans_1d module procedure med_methods_check_for_nans_2d + module procedure med_methods_check_for_nans_3d end interface med_methods_check_for_nans ! used/reused in module @@ -53,8 +55,9 @@ module med_methods_mod public med_methods_FB_getNumflds public med_methods_FB_Field_diagnose public med_methods_FB_GeomPrint - public med_methods_FB_getdata2d public med_methods_FB_getdata1d + public med_methods_FB_getdata2d + public med_methods_FB_getdata3d public med_methods_FB_getmesh public med_methods_FB_check_for_nans @@ -67,6 +70,7 @@ module med_methods_mod public med_methods_Field_getdata1d public med_methods_Field_getdata2d + public med_methods_Field_getdata3d public med_methods_Field_diagnose public med_methods_Field_GeomPrint public med_methods_FieldPtr_compare @@ -113,10 +117,11 @@ subroutine med_methods_FB_init_pointer(StateIn, FBout, flds_scalar_name, name, r integer :: lrank integer :: fieldCount integer :: ungriddedCount - integer :: ungriddedLBound(1) - integer :: ungriddedUBound(1) + integer :: ungriddedLBound(2) + integer :: ungriddedUBound(2) real(R8), pointer :: dataptr1d(:) real(R8), pointer :: dataptr2d(:,:) + real(R8), pointer :: dataptr3d(:,:,:) character(ESMF_MAXSTR), allocatable :: lfieldNameList(:) character(len=*), parameter :: subname='(med_methods_FB_init_pointer)' ! ---------------------------------------------- @@ -161,15 +166,15 @@ subroutine med_methods_FB_init_pointer(StateIn, FBout, flds_scalar_name, name, r call ESMF_FieldGet(lfield, rank=lrank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - if (lrank == 2) then + if (lrank == 3) then ! determine ungridded lower and upper bounds for lfield call ESMF_AttributeGet(lfield, name="UngriddedLBound", convention="NUOPC", & purpose="Instance", itemCount=ungriddedCount, isPresent=isPresent, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - if (ungriddedCount /= 1) then + if (ungriddedCount /= 2) then call shr_log_error(trim(subname)//": ERROR ungriddedCount for "// & - trim(lfieldnamelist(n))//" must be 1 if rank is 2 ", rc=rc) + trim(lfieldnamelist(n))//" must be 2 if rank is 3 ", rc=rc) return end if @@ -181,6 +186,36 @@ subroutine med_methods_FB_init_pointer(StateIn, FBout, flds_scalar_name, name, r purpose="Instance", valueList=ungriddedUBound, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + ! get 3d pointer for field + call ESMF_FieldGet(lfield, farrayptr=dataptr3d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! create new field with two ungridded dimensions + newfield = ESMF_FieldCreate(lmesh, dataptr3d, ESMF_INDEX_DELOCAL, & + meshloc=meshloc, name=lfieldNameList(n), & + ungriddedLbound=ungriddedLbound, ungriddedUbound=ungriddedUbound, gridToFieldMap=(/3/), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + else if (lrank == 2) then + + ! determine ungridded lower and upper bounds for lfield + call ESMF_AttributeGet(lfield, name="UngriddedLBound", convention="NUOPC", & + purpose="Instance", itemCount=ungriddedCount, isPresent=isPresent, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ungriddedCount /= 1) then + call shr_log_error(trim(subname)//": ERROR ungriddedCount for "// & + trim(lfieldnamelist(n))//" must be 1 if rank is 2 ", rc=rc) + return + end if + + ! set ungridded dimensions for field + call ESMF_AttributeGet(lfield, name="UngriddedLBound", convention="NUOPC", & + purpose="Instance", valueList=ungriddedLBound(1:1), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_AttributeGet(lfield, name="UngriddedUBound", convention="NUOPC", & + purpose="Instance", valueList=ungriddedUBound(1:1), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + ! get 2d pointer for field call ESMF_FieldGet(lfield, farrayptr=dataptr2d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return @@ -188,7 +223,7 @@ subroutine med_methods_FB_init_pointer(StateIn, FBout, flds_scalar_name, name, r ! create new field with an ungridded dimension newfield = ESMF_FieldCreate(lmesh, dataptr2d, ESMF_INDEX_DELOCAL, & meshloc=meshloc, name=lfieldNameList(n), & - ungriddedLbound=ungriddedLbound, ungriddedUbound=ungriddedUbound, gridToFieldMap=(/2/), rc=rc) + ungriddedLbound=ungriddedLbound(1:1), ungriddedUbound=ungriddedUbound(1:1), gridToFieldMap=(/2/), rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return else if (lrank == 1) then @@ -203,7 +238,7 @@ subroutine med_methods_FB_init_pointer(StateIn, FBout, flds_scalar_name, name, r if (chkerr(rc,__LINE__,u_FILE_u)) return else - call shr_log_error(trim(subname)//": ERROR only rank1 and rank2 are supported for rank of fields ", rc=rc) + call shr_log_error(trim(subname)//": ERROR only rank1, rank2, and rank3 are supported for rank of fields ", rc=rc) return end if @@ -583,6 +618,7 @@ subroutine med_methods_FB_reset(FB, value, rc) integer :: lrank real(R8), pointer :: fldptr1(:) real(R8), pointer :: fldptr2(:,:) + real(R8), pointer :: fldptr3(:,:,:) character(len=*),parameter :: subname='(med_methods_FB_reset)' ! ---------------------------------------------- @@ -609,7 +645,7 @@ subroutine med_methods_FB_reset(FB, value, rc) do n = 1, fieldCount call ESMF_FieldBundleGet(FB, fieldName=trim(lfieldnamelist(n)), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call med_methods_Field_GetFldPtr(lfield, fldptr1=fldptr1, fldptr2=fldptr2, rank=lrank, rc=rc) + call med_methods_Field_GetFldPtr(lfield, fldptr1=fldptr1, fldptr2=fldptr2, fldptr3=fldptr3, rank=lrank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (lrank == 0) then @@ -618,6 +654,8 @@ subroutine med_methods_FB_reset(FB, value, rc) fldptr1 = lvalue elseif (lrank == 2) then fldptr2 = lvalue + elseif (lrank == 3) then + fldptr3 = lvalue else call shr_log_error(trim(subname)//": ERROR in rank "//trim(lfieldnamelist(n)), & line=__LINE__, file=u_FILE_u, rc=rc) @@ -660,6 +698,7 @@ subroutine med_methods_State_reset(State, value, rc) integer :: lrank real(R8), pointer :: fldptr1(:) real(R8), pointer :: fldptr2(:,:) + real(R8), pointer :: fldptr3(:,:,:) character(len=*),parameter :: subname='(med_methods_State_reset)' ! ---------------------------------------------- @@ -681,7 +720,7 @@ subroutine med_methods_State_reset(State, value, rc) do n = 1, fieldCount call ESMF_StateGet(State, itemName=trim(lfieldnamelist(n)), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call med_methods_Field_GetFldPtr(lfield, fldptr1=fldptr1, fldptr2=fldptr2, rank=lrank, rc=rc) + call med_methods_Field_GetFldPtr(lfield, fldptr1=fldptr1, fldptr2=fldptr2, fldptr3=fldptr3, rank=lrank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (lrank == 0) then ! no local data @@ -689,6 +728,8 @@ subroutine med_methods_State_reset(State, value, rc) fldptr1 = lvalue elseif (lrank == 2) then fldptr2 = lvalue + elseif (lrank == 3) then + fldptr3 = lvalue else call shr_log_error(trim(subname)//": ERROR in rank "//trim(lfieldnamelist(n)), & line=__LINE__, file=u_FILE_u, rc=rc) @@ -720,11 +761,12 @@ subroutine med_methods_FB_average(FB, count, rc) integer , intent(out) :: rc ! local variables - integer :: i,j,n + integer :: i,j,k,n integer :: fieldCount, lrank character(ESMF_MAXSTR) ,pointer :: lfieldnamelist(:) real(R8), pointer :: dataPtr1(:) real(R8), pointer :: dataPtr2(:,:) + real(R8), pointer :: dataPtr3(:,:,:) type(ESMF_Field) :: lfield character(len=*),parameter :: subname='(med_methods_FB_average)' ! ---------------------------------------------- @@ -753,7 +795,7 @@ subroutine med_methods_FB_average(FB, count, rc) do n = 1, fieldCount call ESMF_FieldBundleGet(FB, fieldName=trim(lfieldnamelist(n)), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call med_methods_Field_GetFldPtr(lfield, fldptr1=dataptr1, fldptr2=dataptr2, rank=lrank, rc=rc) + call med_methods_Field_GetFldPtr(lfield, fldptr1=dataptr1, fldptr2=dataptr2, fldptr3=dataptr3, rank=lrank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (lrank == 0) then @@ -768,6 +810,14 @@ subroutine med_methods_FB_average(FB, count, rc) dataptr2(i,j) = dataptr2(i,j) / real(count, R8) enddo enddo + elseif (lrank == 3) then + do k=lbound(dataptr3,3),ubound(dataptr3,3) + do j=lbound(dataptr3,2),ubound(dataptr3,2) + do i=lbound(dataptr3,1),ubound(dataptr3,1) + dataptr3(i,j,k) = dataptr3(i,j,k) / real(count, R8) + enddo + enddo + enddo else call shr_log_error(trim(subname)//": ERROR rank not supported ", rc=rc) return @@ -803,6 +853,7 @@ subroutine med_methods_FB_diagnose(FB, string, rc) character(len=CL) :: lstring real(R8), pointer :: dataPtr1d(:) real(R8), pointer :: dataPtr2d(:,:) + real(R8), pointer :: dataPtr3d(:,:,:) type(ESMF_Field) :: lfield character(len=*), parameter :: subname='(med_methods_FB_diagnose)' ! ---------------------------------------------- @@ -828,7 +879,7 @@ subroutine med_methods_FB_diagnose(FB, string, rc) do n = 1, fieldCount call ESMF_FieldBundleGet(FB, fieldName=trim(lfieldnamelist(n)), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call med_methods_Field_GetFldPtr(lfield, fldptr1=dataptr1d, fldptr2=dataptr2d, rank=lrank, rc=rc) + call med_methods_Field_GetFldPtr(lfield, fldptr1=dataptr1d, fldptr2=dataptr2d, fldptr3=dataptr3d, rank=lrank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (lrank == 0) then @@ -851,6 +902,15 @@ subroutine med_methods_FB_diagnose(FB, string, rc) " no data" endif + elseif (lrank == 3) then + if (size(dataPtr3d) > 0) then + write(msgString,'(A,3g14.7,i8)') trim(subname)//' '//trim(lstring)//': '//trim(lfieldnamelist(n))//' ', & + minval(dataPtr3d), maxval(dataPtr3d), sum(dataPtr3d), size(dataPtr3d) + else + write(msgString,'(A,a)') trim(subname)//' '//trim(lstring)//': '//trim(lfieldnamelist(n)), & + " no data" + endif + else call shr_log_error(trim(subname)//": ERROR rank not supported ", rc=rc) return @@ -1001,6 +1061,7 @@ subroutine med_methods_State_diagnose(State, string, rc) character(len=CS) :: lstring real(R8), pointer :: dataPtr1d(:) real(R8), pointer :: dataPtr2d(:,:) + real(R8), pointer :: dataPtr3d(:,:,:) type(ESMF_Field) :: lfield character(len=*),parameter :: subname='(med_methods_State_diagnose)' ! ---------------------------------------------- @@ -1022,7 +1083,7 @@ subroutine med_methods_State_diagnose(State, string, rc) do n = 1, fieldCount call ESMF_StateGet(State, itemName=trim(lfieldnamelist(n)), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call med_methods_Field_GetFldPtr(lfield, fldptr1=dataptr1d, fldptr2=dataptr2d, rank=lrank, rc=rc) + call med_methods_Field_GetFldPtr(lfield, fldptr1=dataptr1d, fldptr2=dataptr2d, fldptr3=dataptr3d, rank=lrank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (lrank == 0) then @@ -1043,6 +1104,14 @@ subroutine med_methods_State_diagnose(State, string, rc) write(msgString,'(A,a)') trim(subname)//' '//trim(lstring)//': '//trim(lfieldnamelist(n)), & " no data" endif + elseif (lrank == 3) then + if (size(dataPtr3d) > 0) then + write(msgString,'(A,3g14.7,i8)') trim(subname)//' '//trim(lstring)//': '//trim(lfieldnamelist(n)), & + minval(dataPtr3d), maxval(dataPtr3d), sum(dataPtr3d), size(dataPtr3d) + else + write(msgString,'(A,a)') trim(subname)//' '//trim(lstring)//': '//trim(lfieldnamelist(n)), & + " no data" + endif else call shr_log_error(trim(subname)//": ERROR rank not supported ", rc=rc) return @@ -1080,8 +1149,9 @@ subroutine med_methods_FB_Field_diagnose(FB, fieldname, string, rc) character(len=CS) :: lstring real(R8), pointer :: dataPtr1d(:) real(R8), pointer :: dataPtr2d(:,:) + real(R8), pointer :: dataPtr3d(:,:,:) type(ESMF_Field) :: lfield - integer :: ungriddedUBound(1) ! currently the size must equal 1 for rank 2 fields + integer :: fieldrank character(len=*),parameter :: subname='(med_methods_FB_Field_diagnose)' ! ---------------------------------------------- @@ -1097,9 +1167,18 @@ subroutine med_methods_FB_Field_diagnose(FB, fieldname, string, rc) call ESMF_FieldBundleGet(FB, fieldName=fieldname, field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, rc=rc) + call ESMF_FieldGet(lfield, rank=fieldrank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - if (ungriddedUBound(1) > 0) then + if (fieldrank == 3) then + call ESMF_FieldGet(lfield, farrayptr=dataptr3d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (size(dataptr3d) > 0) then + write(msgString,'(A,3g14.7,i8)') trim(subname)//' '//trim(lstring)//': '//trim(fieldname), & + minval(dataPtr3d), maxval(dataPtr3d), sum(dataPtr3d), size(dataPtr3d) + else + write(msgString,'(A,a)') trim(subname)//' '//trim(lstring)//': '//trim(fieldname)," no data" + endif + elseif (fieldrank == 2) then call ESMF_FieldGet(lfield, farrayptr=dataptr2d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (size(dataptr2d) > 0) then @@ -1108,7 +1187,7 @@ subroutine med_methods_FB_Field_diagnose(FB, fieldname, string, rc) else write(msgString,'(A,a)') trim(subname)//' '//trim(lstring)//': '//trim(fieldname)," no data" endif - else + elseif (fieldrank == 1) then call ESMF_FieldGet(lfield, farrayptr=dataptr1d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (size(dataPtr1d) > 0) then @@ -1117,6 +1196,10 @@ subroutine med_methods_FB_Field_diagnose(FB, fieldname, string, rc) else write(msgString,'(A,a)') trim(subname)//' '//trim(lstring)//': '//trim(fieldname)," no data" endif + else + call shr_log_error(subname//": ERROR: unhandled field rank", & + line=__LINE__, file=u_FILE_u, rc=rc) + return end if call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) @@ -1147,6 +1230,7 @@ subroutine med_methods_Field_diagnose(field, fieldname, string, rc) character(len=CS) :: lstring real(R8), pointer :: dataPtr1d(:) real(R8), pointer :: dataPtr2d(:,:) + real(R8), pointer :: dataPtr3d(:,:,:) character(len=*),parameter :: subname='(med_methods_Field_diagnose)' ! ---------------------------------------------- @@ -1183,6 +1267,15 @@ subroutine med_methods_Field_diagnose(field, fieldname, string, rc) else write(msgString,'(A,a)') trim(subname)//' '//trim(lstring)//': '//trim(fieldname)," no data" endif + elseif (lrank == 3) then + call ESMF_FieldGet(field, farrayPtr=dataPtr3d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (size(dataPtr3d) > 0) then + write(msgString,'(A,3g14.7,i8)') trim(subname)//' '//trim(lstring)//': '//trim(fieldname), & + minval(dataPtr3d), maxval(dataPtr3d), sum(dataPtr3d), size(dataPtr3d) + else + write(msgString,'(A,a)') trim(subname)//' '//trim(lstring)//': '//trim(fieldname)," no data" + endif else call shr_log_error(trim(subname)//": ERROR rank not supported ", rc=rc) return @@ -1243,7 +1336,7 @@ subroutine med_methods_FB_accum(FBout, FBin, copy, rc) integer , intent(out) :: rc ! local variables - integer :: i,j,n + integer :: i,j,k,n integer :: fieldCount, lranki, lranko character(ESMF_MAXSTR) ,pointer :: lfieldnamelist(:) logical :: exists @@ -1252,6 +1345,8 @@ subroutine med_methods_FB_accum(FBout, FBin, copy, rc) real(R8), pointer :: dataPtro1(:) real(R8), pointer :: dataPtri2(:,:) real(R8), pointer :: dataPtro2(:,:) + real(R8), pointer :: dataPtri3(:,:,:) + real(R8), pointer :: dataPtro3(:,:,:) type(ESMF_Field) :: lfield character(len=*), parameter :: subname='(med_methods_FB_accum)' ! ---------------------------------------------- @@ -1278,12 +1373,12 @@ subroutine med_methods_FB_accum(FBout, FBin, copy, rc) if (exists) then call ESMF_FieldBundleGet(FBin, fieldName=trim(lfieldnamelist(n)), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call med_methods_Field_GetFldPtr(lfield, fldptr1=dataptri1, fldptr2=dataptri2, rank=lranki, rc=rc) + call med_methods_Field_GetFldPtr(lfield, fldptr1=dataptri1, fldptr2=dataptri2, fldptr3=dataptri3, rank=lranki, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_FieldBundleGet(FBout, fieldName=trim(lfieldnamelist(n)), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call med_methods_Field_GetFldPtr(lfield, fldptr1=dataptro1, fldptr2=dataptro2, rank=lranko, rc=rc) + call med_methods_Field_GetFldPtr(lfield, fldptr1=dataptro1, fldptr2=dataptro2, fldptr3=dataptro3, rank=lranko, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (lranki == 0 .and. lranko == 0) then @@ -1327,6 +1422,31 @@ subroutine med_methods_FB_accum(FBout, FBin, copy, rc) enddo endif + elseif (lranki == 3 .and. lranko == 3) then + + if (.not.med_methods_FieldPtr_Compare(dataPtro3, dataPtri3, subname, rc)) then + call shr_log_error(trim(subname)//": ERROR in dataPtr3 size ", rc=rc) + return + endif + + if (lcopy) then + do k=lbound(dataPtri3,3),ubound(dataPtri3,3) + do j=lbound(dataPtri3,2),ubound(dataPtri3,2) + do i=lbound(dataPtri3,1),ubound(dataPtri3,1) + dataPtro3(i,j,k) = dataPtri3(i,j,k) + enddo + enddo + enddo + else + do k=lbound(dataPtri3,3),ubound(dataPtri3,3) + do j=lbound(dataPtri3,2),ubound(dataPtri3,2) + do i=lbound(dataPtri3,1),ubound(dataPtri3,1) + dataPtro3(i,j,k) = dataPtro3(i,j,k) + dataPtri3(i,j,k) + enddo + enddo + enddo + endif + else write(msgString,'(a,2i8)') trim(subname)//": ranki, ranko = ",lranki,lranko call shr_log_error(trim(subname)//": ERROR ranki ranko not supported "//trim(msgstring)//"\n"//trim(lfieldnamelist(n)), rc=rc) @@ -1395,12 +1515,12 @@ end function med_methods_FB_FldChk !----------------------------------------------------------------------------- - subroutine med_methods_Field_GetFldPtr(field, fldptr1, fldptr2, rank, abort, rc) + subroutine med_methods_Field_GetFldPtr(field, fldptr1, fldptr2, fldptr3, rank, abort, rc) ! ---------------------------------------------- - ! for a field, determine rank and return fldptr1 or fldptr2 + ! for a field, determine rank and return fldptr1, fldptr2, or fldptr3 ! abort is true by default and will abort if fldptr is not yet allocated in field - ! rank returns 0, 1, or 2. 0 means fldptr not allocated and abort=false + ! rank returns 0, 1, 2, or 3. 0 means fldptr not allocated and abort=false ! ---------------------------------------------- use ESMF , only : ESMF_Field,ESMF_Mesh, ESMF_FieldGet, ESMF_MeshGet @@ -1411,6 +1531,7 @@ subroutine med_methods_Field_GetFldPtr(field, fldptr1, fldptr2, rank, abort, rc) type(ESMF_Field) , intent(in) :: field real(R8), pointer , intent(inout), optional :: fldptr1(:) real(R8), pointer , intent(inout), optional :: fldptr2(:,:) + real(R8), pointer , intent(inout), optional :: fldptr3(:,:,:) integer , intent(out) , optional :: rank logical , intent(in) , optional :: abort integer , intent(out) , optional :: rc @@ -1495,6 +1616,15 @@ subroutine med_methods_Field_GetFldPtr(field, fldptr1, fldptr2, rank, abort, rc) call ESMF_FieldGet(field, farrayPtr=fldptr2, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + elseif (lrank == 3) then + if (.not.present(fldptr3)) then + call shr_log_error(trim(subname)//": ERROR missing rank=3 array ", & + line=__LINE__, file=u_FILE_u, rc=rc) + return + endif + call ESMF_FieldGet(field, farrayPtr=fldptr3, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + else call shr_log_error(trim(subname)//": ERROR in rank ", & line=__LINE__, file=u_FILE_u, rc=rc) @@ -1515,7 +1645,7 @@ end subroutine med_methods_Field_GetFldPtr !----------------------------------------------------------------------------- - subroutine med_methods_FB_GetFldPtr(FB, fldname, fldptr1, fldptr2, rank, field, rc) + subroutine med_methods_FB_GetFldPtr(FB, fldname, fldptr1, fldptr2, fldptr3, rank, field, rc) use ESMF , only : ESMF_FieldBundle, ESMF_FieldBundleGet, ESMF_Field @@ -1527,6 +1657,7 @@ subroutine med_methods_FB_GetFldPtr(FB, fldname, fldptr1, fldptr2, rank, field, character(len=*) , intent(in) :: fldname real(R8), pointer , intent(inout), optional :: fldptr1(:) real(R8), pointer , intent(inout), optional :: fldptr2(:,:) + real(R8), pointer , intent(inout), optional :: fldptr3(:,:,:) integer , intent(out), optional :: rank integer , intent(out), optional :: rc type(ESMF_Field) , intent(out), optional :: field @@ -1557,7 +1688,7 @@ subroutine med_methods_FB_GetFldPtr(FB, fldname, fldptr1, fldptr2, rank, field, call ESMF_FieldBundleGet(FB, fieldName=trim(fldname), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call med_methods_Field_GetFldPtr(lfield, fldptr1=fldptr1, fldptr2=fldptr2, rank=lrank, rc=rc) + call med_methods_Field_GetFldPtr(lfield, fldptr1=fldptr1, fldptr2=fldptr2, fldptr3=fldptr3, rank=lrank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (present(rank)) then @@ -1645,6 +1776,43 @@ end function med_methods_FieldPtr_Compare2 !----------------------------------------------------------------------------- + logical function med_methods_FieldPtr_Compare3(fldptr1, fldptr2, cstring, rc) + + ! input/output variables + real(R8), pointer , intent(in) :: fldptr1(:,:,:) + real(R8), pointer , intent(in) :: fldptr2(:,:,:) + character(len=*) , intent(in) :: cstring + integer , intent(out) :: rc + + ! local variables + character(len=*), parameter :: subname='(med_methods_FieldPtr_Compare3)' + ! ---------------------------------------------- + + if (dbug_flag > 10) then + call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) + endif + rc = ESMF_SUCCESS + + med_methods_FieldPtr_Compare3 = .false. + if (lbound(fldptr2,3) /= lbound(fldptr1,3) .or. lbound(fldptr2,2) /= lbound(fldptr1,2) .or. & + lbound(fldptr2,1) /= lbound(fldptr1,1) .or. ubound(fldptr2,3) /= ubound(fldptr1,3) .or. & + ubound(fldptr2,2) /= ubound(fldptr1,2) .or. ubound(fldptr2,1) /= ubound(fldptr1,1)) then + write(msgString,*) trim(subname)//': fldptr2 ',lbound(fldptr2),ubound(fldptr2),': fldptr1 ',lbound(fldptr1),ubound(fldptr1),& + ": ERROR in data size "//trim(cstring) + call shr_log_error(trim(msgString),rc=rc) + return + else + med_methods_FieldPtr_Compare3 = .true. + endif + + if (dbug_flag > 10) then + call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) + endif + + end function med_methods_FieldPtr_Compare3 + + !----------------------------------------------------------------------------- + subroutine med_methods_State_GeomPrint(state, string, rc) use ESMF, only : ESMF_State, ESMF_Field, ESMF_StateGet @@ -1743,6 +1911,7 @@ subroutine med_methods_Field_GeomPrint(field, string, rc) integer :: lrank real(R8), pointer :: dataPtr1(:) real(R8), pointer :: dataPtr2(:,:) + real(R8), pointer :: dataPtr3(:,:,:) type(ESMF_GeomType_Flag) :: geomtype character(len=*),parameter :: subname='(med_methods_Field_GeomPrint)' ! ---------------------------------------------- @@ -1775,7 +1944,7 @@ subroutine med_methods_Field_GeomPrint(field, string, rc) endif call med_methods_Field_GetFldPtr(field, & - fldptr1=dataPtr1, fldptr2=dataPtr2, rank=lrank, abort=.false., rc=rc) + fldptr1=dataPtr1, fldptr2=dataPtr2, fldptr3=dataPtr3, rank=lrank, abort=.false., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (lrank == 0) then @@ -1788,6 +1957,13 @@ subroutine med_methods_Field_GeomPrint(field, string, rc) call ESMF_LogWrite(msgString, ESMF_LOGMSG_INFO) write (msgString,*) trim(subname)//":"//trim(string)//": dataptr bounds dim=2 ",lbound(dataptr2,2),ubound(dataptr2,2) call ESMF_LogWrite(msgString, ESMF_LOGMSG_INFO) + elseif (lrank == 3) then + write (msgString,*) trim(subname)//":"//trim(string)//": dataptr bounds dim=1 ",lbound(dataptr3,1),ubound(dataptr3,1) + call ESMF_LogWrite(msgString, ESMF_LOGMSG_INFO) + write (msgString,*) trim(subname)//":"//trim(string)//": dataptr bounds dim=2 ",lbound(dataptr3,2),ubound(dataptr3,2) + call ESMF_LogWrite(msgString, ESMF_LOGMSG_INFO) + write (msgString,*) trim(subname)//":"//trim(string)//": dataptr bounds dim=3 ",lbound(dataptr3,3),ubound(dataptr3,3) + call ESMF_LogWrite(msgString, ESMF_LOGMSG_INFO) elseif (lrank == 0) then ! means data allocation does not exist yet continue @@ -2309,6 +2485,29 @@ subroutine med_methods_FB_getNumFlds(FB, string, nflds, rc) end subroutine med_methods_FB_getNumFlds + !----------------------------------------------------------------------------- + subroutine med_methods_FB_getdata1d(FB, fieldname, dataptr1d, rc) + + use ESMF, only : ESMF_FieldBundle, ESMF_FieldBundleGet, ESMF_Field, ESMF_FieldGet + + ! input/output variables + type(ESMF_FieldBundle) , intent(in) :: FB + character(len=*) , intent(in) :: fieldname + real(r8) , pointer :: dataptr1d(:) + integer , intent(inout) :: rc + + ! local variables + type(ESMF_Field) :: lfield + ! ---------------------------------------------- + rc = ESMF_SUCCESS + + call ESMF_FieldBundleGet(FB, fieldname=trim(fieldname), field=lfield, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(lfield, farrayptr=dataptr1d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + end subroutine med_methods_FB_getdata1d + !----------------------------------------------------------------------------- subroutine med_methods_FB_getdata2d(FB, fieldname, dataptr2d, rc) @@ -2333,14 +2532,14 @@ subroutine med_methods_FB_getdata2d(FB, fieldname, dataptr2d, rc) end subroutine med_methods_FB_getdata2d !----------------------------------------------------------------------------- - subroutine med_methods_FB_getdata1d(FB, fieldname, dataptr1d, rc) + subroutine med_methods_FB_getdata3d(FB, fieldname, dataptr3d, rc) use ESMF, only : ESMF_FieldBundle, ESMF_FieldBundleGet, ESMF_Field, ESMF_FieldGet ! input/output variables type(ESMF_FieldBundle) , intent(in) :: FB character(len=*) , intent(in) :: fieldname - real(r8) , pointer :: dataptr1d(:) + real(r8) , pointer :: dataptr3d(:,:,:) integer , intent(inout) :: rc ! local variables @@ -2350,10 +2549,27 @@ subroutine med_methods_FB_getdata1d(FB, fieldname, dataptr1d, rc) call ESMF_FieldBundleGet(FB, fieldname=trim(fieldname), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_FieldGet(lfield, farrayptr=dataptr1d, rc=rc) + call ESMF_FieldGet(lfield, farrayptr=dataptr3d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - end subroutine med_methods_FB_getdata1d + end subroutine med_methods_FB_getdata3d + + !----------------------------------------------------------------------------- + subroutine med_methods_Field_getdata1d(field, dataptr1d, rc) + + use ESMF, only : ESMF_Field, ESMF_FieldGet + + ! input/output variables + type(ESMF_Field) , intent(in) :: field + real(r8) , pointer :: dataptr1d(:) + integer , intent(inout) :: rc + ! ---------------------------------------------- + rc = ESMF_SUCCESS + + call ESMF_FieldGet(field, farrayptr=dataptr1d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + end subroutine med_methods_Field_getdata1d !----------------------------------------------------------------------------- subroutine med_methods_Field_getdata2d(field, dataptr2d, rc) @@ -2373,21 +2589,21 @@ subroutine med_methods_Field_getdata2d(field, dataptr2d, rc) end subroutine med_methods_Field_getdata2d !----------------------------------------------------------------------------- - subroutine med_methods_Field_getdata1d(field, dataptr1d, rc) + subroutine med_methods_Field_getdata3d(field, dataptr3d, rc) use ESMF, only : ESMF_Field, ESMF_FieldGet ! input/output variables type(ESMF_Field) , intent(in) :: field - real(r8) , pointer :: dataptr1d(:) + real(r8) , pointer :: dataptr3d(:,:,:) integer , intent(inout) :: rc ! ---------------------------------------------- rc = ESMF_SUCCESS - call ESMF_FieldGet(field, farrayptr=dataptr1d, rc=rc) + call ESMF_FieldGet(field, farrayptr=dataptr3d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - end subroutine med_methods_Field_getdata1d + end subroutine med_methods_Field_getdata3d !----------------------------------------------------------------------------- subroutine med_methods_FB_getmesh(FB, mesh, rc) @@ -2433,6 +2649,7 @@ subroutine med_methods_FB_check_for_nans(FB, maintask, logunit, rc) character(len=CL) :: fieldname real(r8) , pointer :: dataptr1d(:) real(r8) , pointer :: dataptr2d(:,:) + real(r8) , pointer :: dataptr3d(:,:,:) integer :: nancount character(len=CS) :: nancount_char character(len=CL) :: msg_error @@ -2459,10 +2676,18 @@ subroutine med_methods_FB_check_for_nans(FB, maintask, logunit, rc) call ESMF_FieldGet(field, farrayPtr=dataptr1d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call med_methods_check_for_nans(dataptr1d, nancount) - else + elseif (fieldrank == 2) then call ESMF_FieldGet(field, farrayPtr=dataptr2d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call med_methods_check_for_nans(dataptr2d, nancount) + elseif (fieldrank == 3) then + call ESMF_FieldGet(field, farrayPtr=dataptr3d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call med_methods_check_for_nans(dataptr3d, nancount) + else + call shr_log_error(subname//": ERROR: unhandled field rank", & + line=__LINE__, file=u_FILE_u, rc=rc) + return end if if (nancount > 0) then write(nancount_char, '(i0)') nancount @@ -2485,11 +2710,11 @@ subroutine med_methods_check_for_nans_1d(dataptr, nancount) real(r8) , intent(in) :: dataptr(:) integer , intent(out) :: nancount ! local variables - integer :: n + integer :: i nancount = 0 - do n = 1,size(dataptr) - if (shr_infnan_isnan(dataptr(n))) then + do i = 1,size(dataptr) + if (shr_infnan_isnan(dataptr(i))) then nancount = nancount + 1 end if end do @@ -2501,16 +2726,36 @@ subroutine med_methods_check_for_nans_2d(dataptr, nancount) real(r8) , intent(in) :: dataptr(:,:) integer , intent(out) :: nancount ! local variables - integer :: n,k + integer :: i,j nancount = 0 - do k = 1,size(dataptr, dim=1) - do n = 1,size(dataptr, dim=2) - if (shr_infnan_isnan(dataptr(k,n))) then + do j = 1,size(dataptr, dim=2) + do i = 1,size(dataptr, dim=1) + if (shr_infnan_isnan(dataptr(i,j))) then nancount = nancount + 1 end if end do end do end subroutine med_methods_check_for_nans_2d + subroutine med_methods_check_for_nans_3d(dataptr, nancount) + use shr_infnan_mod, only: shr_infnan_isnan + ! input/output variables + real(r8) , intent(in) :: dataptr(:,:,:) + integer , intent(out) :: nancount + ! local variables + integer :: i,j,k + + nancount = 0 + do k = 1,size(dataptr, dim=3) + do j = 1,size(dataptr, dim=2) + do i = 1,size(dataptr, dim=1) + if (shr_infnan_isnan(dataptr(i,j,k))) then + nancount = nancount + 1 + end if + end do + end do + end do + end subroutine med_methods_check_for_nans_3d + end module med_methods_mod diff --git a/mediator/med_phases_prep_glc_mod.F90 b/mediator/med_phases_prep_glc_mod.F90 index 5007801c..f3cbad49 100644 --- a/mediator/med_phases_prep_glc_mod.F90 +++ b/mediator/med_phases_prep_glc_mod.F90 @@ -752,7 +752,6 @@ subroutine med_phases_prep_glc_map_lnd2glc(gcomp, rc) type(InternalState) :: is_local real(r8), pointer :: topolnd_g_ec(:,:) ! topo in elevation classes real(r8), pointer :: topoglc_g(:) ! ice topographic height on the glc grid extracted from glc import - real(r8), pointer :: data_ice_covered_g(:) ! data for ice-covered regions on the GLC grid real(r8), pointer :: ice_covered_g(:) ! if points on the glc grid is ice-covered (1) or ice-free (0) integer , pointer :: elevclass_g(:) ! elevation classes glc grid real(r8), pointer :: dataexp_g(:) ! pointer into @@ -760,6 +759,10 @@ subroutine med_phases_prep_glc_map_lnd2glc(gcomp, rc) real(r8) :: elev_l, elev_u ! lower and upper elevations in interpolation range real(r8) :: d_elev ! elev_u - elev_l integer :: nfld, ec + integer, allocatable :: index_lower(:) ! lower EC index for vertical interpolation, per glc gridcell + integer, allocatable :: index_upper(:) ! upper EC index for vertical interpolation, per glc gridcell + real(r8), allocatable :: weight_lower(:) ! weight for lower EC, per glc gridcell + real(r8), allocatable :: weight_upper(:) ! weight for upper EC, per glc gridcell integer :: n,lsize_g,ns type(ESMF_Field) :: field_lfrac_l integer :: fieldCount @@ -866,6 +869,60 @@ subroutine med_phases_prep_glc_map_lnd2glc(gcomp, rc) call fldbun_getdata2d(toglc_frlnd(ns)%FBlndAccum2glc_g, 'Sl_topo_elev', topolnd_g_ec, rc=rc) if (chkErr(rc,__LINE__,u_FILE_u)) return + ! ------------------------------------------------------------------------ + ! Pre-compute vertical interpolation indices and weights for each glc gridcell. + ! These depend only on topography (not on the field being interpolated), so we + ! compute them once here and reuse them for each field below. + ! ------------------------------------------------------------------------ + + allocate(index_lower(lsize_g)) + allocate(index_upper(lsize_g)) + allocate(weight_lower(lsize_g)) + allocate(weight_upper(lsize_g)) + + do n = 1, lsize_g + if (topoglc_g(n) < topolnd_g_ec(2,n)) then + ! lower than lowest mean EC elevation: use lowest EC value + ! (note that index 1 is bare land, so index 2 is lowest EC) + index_lower(n) = 2 + index_upper(n) = 2 + weight_lower(n) = 0._r8 + weight_upper(n) = 1._r8 + else if (topoglc_g(n) >= topolnd_g_ec(ungriddedCount,n)) then + ! higher than highest mean EC elevation: use highest EC value + index_lower(n) = ungriddedCount + index_upper(n) = ungriddedCount + weight_lower(n) = 1._r8 + weight_upper(n) = 0._r8 + else + ! find bounding ECs and linearly interpolate + do ec = 3, ungriddedCount + if (topoglc_g(n) < topolnd_g_ec(ec,n)) then + elev_l = topolnd_g_ec(ec-1,n) + elev_u = topolnd_g_ec(ec ,n) + d_elev = elev_u - elev_l + index_lower(n) = ec-1 + index_upper(n) = ec + if (d_elev <= 0._r8) then + ! This shouldn't happen, but handle it in case it does. In this case, + ! let's arbitrarily use the mean of the two elevation classes, rather + ! than the weighted mean. + write(logunit,*) subname//' WARNING: topo diff between elevation classes <= 0' + write(logunit,*) 'n, ec, elev_l, elev_u = ', n, ec, elev_l, elev_u + write(logunit,*) 'Simply using mean of the two elevation classes,' + write(logunit,*) 'rather than the weighted mean.' + weight_lower(n) = 0.5_r8 + weight_upper(n) = 0.5_r8 + else + weight_lower(n) = (elev_u - topoglc_g(n)) / d_elev + weight_upper(n) = (topoglc_g(n) - elev_l) / d_elev + end if + exit + end if + end do + end if + end do + ! ------------------------------------------------------------------------ ! Loop over fields in export field bundle to glc for ice sheet ns and ! perform vertical interpolation of data onto ice sheet topography @@ -877,7 +934,6 @@ subroutine med_phases_prep_glc_map_lnd2glc(gcomp, rc) ! current glint implementation, which sets acab and artm to 0 over ocean (although ! notes that this could lead to a loss of conservation). Figure out how to handle this case. - allocate(data_ice_covered_g(lsize_g)) do nfld = 1, size(fldnames_to_glc) ! Get a pointer to the land data in multiple elevation classes on the glc grid @@ -888,58 +944,16 @@ subroutine med_phases_prep_glc_map_lnd2glc(gcomp, rc) call fldbun_getdata1d(is_local%wrap%FBExp(compglc(ns)), fldnames_to_glc(nfld), dataexp_g, rc) if (chkErr(rc,__LINE__,u_FILE_u)) return - ! First set data_ice_covered_g to bare land everywehre - data_ice_covered_g(:) = 0._r8 - - ! Loop over land points and overwrite with valid values + ! Apply pre-computed vertical interpolation indices and weights do n = 1, lsize_g - - ! For each ice sheet point, find bounding EC values... - if (topoglc_g(n) < topolnd_g_ec(2,n)) then - - ! lower than lowest mean EC elevation value - data_ice_covered_g(n) = dataptr2d(2,n) - - else if (topoglc_g(n) >= topolnd_g_ec(ungriddedCount, n)) then - - ! higher than highest mean EC elevation value - data_ice_covered_g(n) = dataptr2d(ungriddedCount,n) - - else - - ! do linear interpolation of data in the vertical - do ec = 3, ungriddedCount - if (topoglc_g(n) < topolnd_g_EC(ec,n)) then - elev_l = topolnd_g_EC(ec-1,n) - elev_u = topolnd_g_EC(ec ,n) - d_elev = elev_u - elev_l - if (d_elev <= 0) then - ! This shouldn't happen, but handle it in case it does. In this case, - ! let's arbitrarily use the mean of the two elevation classes, rather - ! than the weighted mean. - write(logunit,*) subname//' WARNING: topo diff between elevation classes <= 0' - write(logunit,*) 'n, ec, elev_l, elev_u = ', n, ec, elev_l, elev_u - write(logunit,*) 'Simply using mean of the two elevation classes,' - write(logunit,*) 'rather than the weighted mean.' - data_ice_covered_g(n) = dataptr2d(ec-1,n) * 0.5_r8 & - + dataptr2d(ec ,n) * 0.5_r8 - else - data_ice_covered_g(n) = dataptr2d(ec-1,n) * ((elev_u - topoglc_g(n)) / d_elev) & - + dataptr2d(ec ,n) * ((topoglc_g(n) - elev_l) / d_elev) - end if - exit - end if - end do - end if ! topoglc_g(n) - if (elevclass_g(n) /= 0) then - ! ice-covered cells have interpolated values - dataexp_g(n) = data_ice_covered_g(n) + ! ice-covered cells: vertically interpolate between bounding ECs + dataexp_g(n) = dataptr2d(index_lower(n), n) * weight_lower(n) & + + dataptr2d(index_upper(n), n) * weight_upper(n) else - ! non ice-covered cells have bare land value - dataexp_g(n) = dataptr2d(1,n) + ! non ice-covered cells: use bare land value (EC index 1) + dataexp_g(n) = dataptr2d(1, n) end if - end do ! end of loop over land points end do ! end loop over fields (nflds) @@ -961,7 +975,10 @@ subroutine med_phases_prep_glc_map_lnd2glc(gcomp, rc) ! clean up memory that is ice sheet dependent deallocate(elevclass_g) - deallocate(data_ice_covered_g) + deallocate(index_lower) + deallocate(index_upper) + deallocate(weight_lower) + deallocate(weight_upper) end do ! end of loop over ice sheets